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Introductlon 

This  paper  discusses  adaptive  meshes  that  are  generated  by  the  use  of 
variational  principles.  These  adaptive  meshes  are  used  in  conjunction  with 
finite  difference  methods  to  continuously  rezone  calculations.  The  fixed  number 
of  mesh  vertices  are  moved  in  a  manner  which  enhances  the  numerical  accuracy  of 
these  finite  difference  schemes  when  solving  singular  problems. 

Adaptive  meshes  are  the  subject  of  a  rather  large  and  varied  research 
effort.  However,  a  definition  can  be  developed  that  encompasses  most  of  the 
areas  of  research.  An  adaptive  computation  mesh  is  one  that  is  generated  in 
response  to  the  structure  of  a  problem  so  as  to  enhance  the  numerical  accuracy  of 
the  associated  differencing  scheme.  It  is  the  method  of  generating  the  mesh  that 
gives  this  area  its  diversity. 

Several  different  methods  of  generating  an  adaptive  mesh  will  be  discussed. 
The  methods  can  be  categorized  in  many  ways,  but  three  categories  are  chosen 
here.   They  are  fixed  meshes,  selectively  refined  meshes  and  moving  meshes. 

A  fixed  mesh  is  one  that  has  the  points  of  the  computation  nesh  fixed  durinp 
the  evolution  or  solution  of  a  problem.  For  example,  in  viscous  flows  boundary 
layers  are  present  at  no-slip  boundaries  and  their  thickness  can  he  predicted 
within  an  order  of  magnitude.  Thompson  et  al  [2A]  and  Sorenson  [17A]  use  this 
knowledge  to  generate  body  fitted  grids  that  place  many  grid  points  close  to  the 
no-slip  boundaries.  They  solve  potential  equations  with  source  terms  to  generate 
a  smooth  mesh  in  the  Interior  while  compressing  the  grid  near  the  boundary.  Sone 
theoretical  work  has  been  done  in  this  area  by  Browning  et  al.  [18A]  and  by 
Ciment  [19A].  Browning  discussed  stability  problems  for  nonuniform  grids  in  one 
dimension  and  Ciment  has  derived  stability  criteria  in  both  one  and  two 
dimensions. 
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A  selectively  refined  mesh  technique  is  a  method  that  adds  or  substracts 
computation  points  in  response  to  how  a  problem  evolves.  The  two-point  boundary 
value  problem  has  been  an  area  of  much  activity  for  selective  refinement. 
Babuska  and  Rheinboldt  [20A]  have  given  a  survey  of  the  literature  and  an 
analysis  of  optimal  finite  element  meshes.  Essentially  most  of  these  methods 
solve  the  boundary  value  problem  more  than  once.  At  each  Iteration  an  estimate 
is  made  of  the  error  of  the  solution.  Then,  points  are  added  or  subtracted  to 
reduce  the  error.  Finally  a  new  solution  is  found  and  the  process  is  repeated 
until  certain  error  criteria  are  met. 

Time  dependent  problems  in  one  dimension  can  use  many  of  the  techniques 
developed  for  the  two-point  boundary  value  problem.  Chong  [21A]  has  developed  a 
refinement  technique  that  estimates  where  the  mesh  should  be  refined  every  few 
cycles  in  the  calculation.  Areas  that  surround  the  refined  mesh  are  also  refined 
to  give  a  safety  margin  until  the  next  refinement  is  done. 

Gropp  [22A]  and  Lam  and  Simpson  [23A]  have  done  selectively  refined  mesh 
calculations  In  two  dimensions.  Gropp  calculated  the  solution  of  the  2-D 
inviscid  time  dependent  Burger  Equation  while  Lam  and  Simpson  solved  a  stearty 
state  advection-dif fusion  equation.  McBryan's  [24A]  front  tracking  algorithm 
also  falls  into  the  category  of  selective  refinement.  This  algorithm  not  only 
refines  a  mesh,  but  aligns  the  mesh  parallel  to  discontinuities  in  the  solution 
for  greater  resolution. 

The  last  class  of  adaptive  meshes  are  the  moving  meshes.  Moving  meshes  are 
meshes  that  enhance  resolution  by  bunching  a  fixed  number  of  points  at  the 
appropriate  places.  Gellnas  and  Miller  [25A]  have  devised  a  finite  element 
formulation  that  Includes  node  locations  as  part  of  the  variational  formulation. 
That  Is  the  nodes  move  and  the  solution  is  altered  simultaneously  to  minimize  the 
error   functional.    Brackbill   and   Saltzman   [9A]   and   Yanenko   [8A]   have 
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Lagrangean  hydrodynamics  calculations  have  been  around  for  a  very  long  time 
[IB].  In  one  dimension  Lagrangean  calculations  have  often  been  quite  accurate. 
However,  many  two-dimensional  flows  tend  to  tangle  a  Lagrangean  mesh  hopelessls . 
The  laboratory  or  Eulerian  frame  equations  circumvent  mesh  tangling  but  introduce 
new  problems  in  accurately  convecting  mass,  momentum,  and  energy  [5A,6A].  An 
alternative  way  of  solving  the  problem  of  mesh  tangling  is  to  rezone  a  Lagrangean 
mesh  only  as  needed.  Hirt  et  al.  [7A]  developed  a  method  that  allows  the  mesh 
to  move  in  an  arbitrary  manner.  This  technique  was  then  used  by  Brackbill  and 
Pracht  [26A]  to  continuously  rezone  the  time  dependent  problem  using  Winslow's 
[27A]  mesh  generator.   This  is  the  other  starting  point  of  this  paper. 

By  recongnizing  that  the  elliptic  systems  used  for  generating  body  fitted 
meshes  can  be  derived  from  a  simple  variational  principle  [28A]  that  measures 
mesh  smoothness,  other  variational  principles  can  be  added  to  make  the  mesh 
adaptive  [8A,9A].  Since  an  adaptive  mesh  is  a  moving  mesh,  Hirt's  formulation  is 
again  a  convenient  device  for  allowing  the  mesh  to  move. 

The  scheme  of  this  thesis  will  be  to  describe  adaptive  mesh  techniques  in 
one  through  three  spatial  dimensions.  The  first  chapter  describes  a  one- 
dimensional  mesh,  gives  some  asymptotic  results  with  a  simple  numerical 
experiment . 

Chapters  two  and  three  will  extend  the  reference  [9A]  by  presenting  both 
alternative  variational  principles  and  extensions  in  higher  dimensions  of 
variational  principles  that  have  been  already  introduced.  In  particular,  chapter 
two  will  deal  with  two-dimensional  meshes  while  chapter  three  will  deal  with 
three-dimensional  meshes. 

The  last  chapter  will  discuss  some  practical  applications  of  adaptive  meshes 
to  two-dimensional  problems.  The  first  problem  to  be  presented  is  thr  classical 
two  fluid  gravitational  instability  (Rayleigh-Taylor)  [5B].   The  second  problem 
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is  the  wind  tunnel  with  a  step  studied  by  Emery,  van  Leer,  Woodward,  and  Coella 
[10A1,11A,12A]. 
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I .   Mesh  Generation  and  Adaptive  Mesh  Techniques  In  1-D 
Computation  Meshes  as  Mappings 

When  solving  a  1-D  partial  differential  equation  by  a  finite  difference 
technique  an  approximate  solution  is  calculated  at  a  set  of  well  ordered  points 

x^,  X2,  .  ••  ,  Xj^  A  xj^  <  Xj    1  <  j   .  (I.l) 

The  subscripts  are  used  by  the  computation  device  explicitly  to  keep  track  of  the 
points  on  the  mesh.  In  addition,  the  subscripts  also  suggest  a  mapping  from  a 
subset  of  the  integers  Z  to  a  subset  of  the  real  line  t. 


Xj:  Z  -^  *   .  (1.2) 


It  will  be  convenient  however  to  think  of  the  mapping  as  one  that  maps  a 
subset  of  the  real  line  onto  another  subset  of  the  line.   Using  I.l  as  an  example 

X-  [1,N]  *  [xi,  x^]   .  (1.3) 

The  domain  of  the  mapping  will  be  called  the  parameter  space  and  the  range 
of  the  mapping  will  be  called  a  physical  space.  Elements  of  the  parameter  space 
will  have  a  variable  name  ^  and  elements  of  the  physical  space  will  have  variable 
name  x.  Therefore  the  mapping  x-  becomes  x(C).  We  will  also  assume  that  x(0  is 
1-1  and  onto  so  that  5(x)  exists. 

Mesh  generation  using  a  variational  formulation 

Suppose  there  is  a  function  w(x)  >  w^  >  0  on  (X|,Xj^)  associated  with  a  ^PF . 
Further,   assume  w   is   large  where  a   large  number  of  points   is  needed   to 
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approxlmate  the  solution  of  the  associated  partial  differential  equation  to  sone 
fixed  error.   Where  w  is  smaller  fewer  points  are  needed  to  get  the  same  error. 

If  given  N  points  to  do  the  finite  difference  computation,  and  w,  where 
should  the  points  x-  be  placed  as  to  minimize  the  error?  To  answer  this  question 
consider 

'Si 

I^(x)  =  /   w(x)x^  dx   .  (I. A) 

'^l  -.     ..      

Here  x  is  an  element  of  C  [1,N]  (continuously  dif f erentiable  functions  on  the 
closed  interval  [1,N]).  Further  restrict  x  to  be  a  strictly  monotone  increasing 
function  (Xf  >  0)  and  x(l)  =  x^ ,  x(N)  =  Xj^. 

To  first  order  Xr  can  be  considered  a  measure  of  the  mesh  spacing  by  the 
mean  value  theorem 


* 

dx(x   )  ,    *   /  /T   CN 


Then  to  first  order,  to  make  Iy(x)  small,  where  w  is  large  Xr  should  be  small. 
When  w  is  small  Xr  need  not  be. 

The  functional  Iy(x)  then  suggests  a  variational  formulation  as  a  method  of 
generating  a  computation  mesh.  That  is  minimize  (I. A)  with  respect  to  the  given 
class  of  functions  and  then  calculate  x(^ )  where  ^  =  1,  ...,  N  to  get  the  mesh. 
To  Fi:m''^i/'e  T   write  (I. A)  in  a  form  where  C  is  dependent  on  x. 

\(K)    =  /    w(x)a^)  ^  dx   .  (T.6) 

^1 

Then  the  Euler  equation  is 
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t^-^aC)  I<w(.)«,)->1  -  0  .  (1.7, 


where 


CCx,)  =  1,  axx,)  =  N 


Since  w(£  )   does  not  explicitly  depend  on  E,    then 


w(C^)-2  =  wx^2  =  ,^ 


or  (T.8) 


w^/\  =  ci 


Separating  variables  and  using  the  boundary  conditions  in  (1.7)  to  deter?iine  the 
two  constants  of  integration 


X 

,1/2, 


/  w^'^is)    ds 
4  =  1  +  (N  -  1)  — .  (1.9) 


Then  (1.9)  can  be  solved  for  x.  for  i  =  1,  ...,  N.  Several  things  should  he 
noted.  First,  from  (1.8)  the  solution  gives  an  equal  error  solution.  That  is 
the  error  would  be  equally  distributed  over  the  mesh.  Using  w  >  w  >  0  and  (1.9) 
we  see  that  x(C  )  is  strictly  monotone  increasing  and  satisfies  the  boundary 
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condition  In  (1.7)  so  we  have  a  solution  to  the  Euler  equation  If  w  Is  smooth 
enough. 

If  the  variation  of  w  is  large  the  mesh  spacing  will  vary  greatly.  Often  it 
is  necessary  to  avoid  this  great  variation  in  mesh  spacing  because  of  a  loss  of 
accuracy  by  the  differencing.  Turning  to  another  functional  we  can  measure  the 
smoothness  of  a  mesh  with 

igCx)  =  /  (C^)2  dx  .  (I. in) 

Solving  the  associated  Euler  equation  of  (I. 10) 


(X  -  Xi ) 

C  =  1  +  (K  -  1) .  (1. 11) 

(X.,  -  Xn) 


This  solution  has  no  variation  in  mesh  spacing  and  really  is  the  smoothest 
mapping  that  could  be  done.  Then,  to  control  the  mesh  spacing  variation  a  linear 
combination  of  the  functionals  is  proposed: 

l(x)  =  (1  -  t)  /   (C^)^  dx  +  t  /   w(x)  x^  dx   .  (1.12) 

x^  ^1 

By  varying  t  between  0  and  1  the  smoothness  of  the  mapping  should  decrease  toward 
the  solution  (1.7). 

Since  the  Euler  equation  will  have  to  be  solved  numerically  at  the  points 
C  =1,  .,.,  N  it  is  best  to  write  (1.12)  using  I,  as  the  independent  variable. 
Since  ^  =  1,2,  ,..,  N  is  uniform  in  spacing  accurate  centered  differences  can  be 
used. 

N  N 

I(x)  =  (1  -  t)  /  (x^)"^  dC  +  t  /  wx^2  d^  (1.13) 

1  1 
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Writlng  out  the  Euler  equation 


^-ia^K'-'H"'^--^'!-" 


(I.IM 


or 


\^    + 


tw 


x'^ 


=  0 


2[(1  -  t)  +  x^^wt] 


(1.15) 


with 


x(l)  =  x^    x(N)  =  xx. 

This  equation  is  nonlinear  and  will  not  yield  a  solution  by  separation  of 
variables.   However  a  great  deal  can  be  said  about  the  solution. 

In  Appendix  A  a  proof  of  the  existence  of  (1.15)  is  given.  In  the  course  of 
the  proof  many  properties  of  the  solution  of  (1.15)  are  illuninated  bv  thp 
estimates.   However,  for  continuity,  we  state  only  the  results  here: 

Theorem  1:   Suppose  w(x)  is  continuously  dif f erentiable  and  bounded  away  fron 
zero: 


w(x)  >  w^  >  0 


Then  the  solution  to  (1.15)  exists  and  is  unique  for  T  e  [0,1 
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Model  Problem  .,  . 

In  this  section  we  will  examine  numerical  approximations  of  a  traveling 
singular  wave  in  one  dimension.  We  will  show  that  an  adaptive  mesh  enhances 
resolution  of  the  singular  region  in  a  predictable  manner.  The  equation  that 
describes  the  wave  is: 


U^  +  CUj^  =  <Uy.^  (C,K   >  0) 


u(l,t)  =  0   u(0,t)  =  1    u(x,0)  =  f(x) 


(1.16) 


Qualitatively  the  solution  will  evolve  in  a  manner  that  moves  the  initial 
data  f(x)  to  the  right  at  a  velocity  c.  In  addition  to  the  translation  the 
sharpness  of  the  initial  data  will  diffuse  in  time.  As  t  *  «>  a  boundary  layer  is 
formed  at  x  ~  1  with  a  thickness  proportional  to  <  . 

With  the  above  in  raind,  the  regime  of  interest  of  (1.16)  will  be  when  k  <<  1 
and  c  =  0(1).  In  this  regime  an  adaptive  mesh  technique  will  be  able  to 
selectively  resolve  a  boundary  layer  without  having  to  waste  points  elsewhere. 

The  solution  of  (1.16)  can  be  written  as  the  sum  of  a  transient  solution 
Uy(x,t)  and  a  steady  state  solution  Ug(x). 

u(x,t)  =  Ug(x)  +  uj(x,t)   .  (1.17) 

Setting  Uj.  =  0  in  (1.16) 

f  1  _  eC(x-l)/'C] 

u_(x)  =  Li -, L      .  (1.18) 

s       (1  -  e-'^/^j 
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For  c/ic  >>  1 ,  a  boundary  layer  near  x  =  1  is  evident  in  Fig.  1.   The  boundary 
layer  1/2  thickness  is  defined  as  the  distance  1  -  x  where  Ug(x)  is  equal  to  1/2. 
6  is  approximately 


6  =  0.69  1   .  .  (1.19) 


The  transient  solution  satisfies  the  PDE  in  (1.16)  but  now  has  different  boundary 
conditions  and  initial  data 

u^d.t)  =  u^(0,t)  =  0   u^(x,0)  =  f(x)  -  Ug(x)   .  (1.20) 

To  solve  the  transient  solution  the  method  of  separation  of  variables  is 
used.   Assume 

u^  =  X(x)T(t)  (1.21) 

and  is  the  sum  of  a  countable  number  of  exponentially  decaying  harmonics 


T„   <X„"   cX„' 
Tn    \  >''n 


where  (*)  =  d/dt  and  (')  =  d/dx.   The  solution  to  the  tiire  (t)  dependent  part  of 
the  equation  is 

-n't 
Tn  =  e   "  (1.23) 

For  the  space  part  of  the  equation,  the  differential  equation  is 
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2 

c       "n 
K"   -  -  X„'  +  J—   X„  =  0 

n    1^   n     <    n 

(I.2A) 
\(0)  =  0   X^(l)  =  0   . 

The  equation  can  be  transformed  to  a  Sturm-Liouville  type  equation  by  ciultirlyinj 
through  by 

r(x)  =  e""^'^/^  (1.25) 

yielding, 


d 
dx  ^ "      dx 


From   general   theory   [12B]   there   exists   a   sequence   of   eigenvalues   and 
corresponding  eigenfunctions  orthogonal  with  respect  to  the  inner  product 


6ij  =  (X^.Xj)  =  /  ^ X^Xj  dx   .  (1.26) 

o 


The   eigenvalues    are 


a^2   =  ^n^Ti^  +  c^/Ak        n   =    1,    2,    ...  (1.27) 


and  the  corresponding  eigenfunctions  are 

\^   =  /2<  e^/^  sin(mx)   n  =  1,  2,  ...   .  (1.28) 
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Using  (1.21)  thru  (1.28)  we  find 

2  2 

e 
n=l 


"   -(<nS^+(cB/4))t     ,/,  „^/, 
u(x,t)  =  I   e  (2<)l/2e6x/2 

(1.29) 


r^  e-^^/2       , 

sln(mTx)  /   g(0  ^ (2»c)^'^  slndrrC)  d?  +  u^(x) 

0 


where 


=  c/<  and  g(5)  =  f(5)  -  Ug(0 


Unfortunately  regimes  of  interest  are  when  B  >>  1.  In  this  regime  the  solution 
is  badly  conditioned  making  (1.29)  useless  as  a  method  for  computation  of  a 
solution.  However,  the  solution  can  be  resummed.  First  the  transient  part  of 
the  solution  (1.29)  may  be  expressed  in  terms  of  the  classical  theta  function 
(See  Abramowitz  and  Stegun  [9B]) 

63(2, q)  =  5;    q"  e2niz  =  1  +  2  ^    q"   cos(2nz)  (1.30) 

n=-"  n=l 


u(x.t)  =  /   g(0  ee/2[(x-^)-ct/2]  i/2[e3(.  ^JL^    ,    e"-'^) 


(1.31) 


-e3[.  1^4i2  .  e--'^)]  6, 
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By  using  the  Poisson  Summation  Formula  (See  Courant-Hilbert  Volume  1,  p.  1  k-11 
[lOB]) 


e(z,e-^^^^)   =     I        e^^^"^'^   e2"iz  =  _i_     [        ^-(0-2/71  )2/Kt      ^ 
n=-<»>  /icTi  t   n=-«' 


Using  the  identity  (1.32)  in  (1.31)  the  solution  may  be  written 


u(x,t)  =  L_./'  e-Kx-ct)-^]2/4ctg(^) 

(AKTIt)^/^  0 


.    [    I        e[v2-v(x-^)]/<t   _  g-[v2-v(x4i)+x4]/<t]    j^ 


(1.32) 


(1.33) 


In  this  form  the  solution  can  be  computed  accurately  using  numerical  quadrature. 
This  is  facilitated  by  the  fact  that  for  <  >>  1  only  a  small  number  of  terms  need 
be  computed  fron  the  series  in  (1.33).  Further,  the  first  term  in  (1.33) 
concentrates  most  of  its  weight  at  the  point  x  -  ct  =  ^ . 

Differentiating  (1.33)  will  yield  a  stable  method  for  computing  derivatives 
of  the  solution.   In  Fig.  2  there  are  some  sample  calculations  using  (1.33)  with 


1     ^  .  X  -  0.5  >, 

tanh  [  — ) 

f(C)  =  i -4 +0.5   .  (1.34) 

tanh  (  M  ) 


There  are  several  important  features  of  the  solution  to  notice.  First  It 
should  be  noted  that  for  Kt  «  1,  the  solution  and  its  derivatives  are  behaving 
like  the  one-dimensional  wave  operator  solution  of 
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u^  +  cu^  =  0   .  (1.35) 


However,  as  k:  t  gets  larger  the  solution  of  (1.16)  exhibits  a  degradation  in 
sharpness  compared  to  the  earlier  time. 

Another  important  feature  to  notice  is  that  all  the  derivatives  become  large 
around  the  advancing  front  of  the  solution  and  are  small  elsewhere.  It  is  this 
feature  of  the  solution  that  makes  an  adaptive  mesh  technique  advantageous. 

A  Finite  Difference  Scheme  for  the  PDE  (1.16) 

Many  numerical  calculations  have  been  successful,  in  part,  because  they 
incorporate  properties  of  a  differential  equation  into  the  approximation.  The 
conservation  laws  for  compressible  inviscid  flows,  for  example,  have  been 
successfully  differenced  in  conservation  form  [IB].  The  one  outstanding  property 
of  (1.16)  is  the  maximum  principle.  That  is,  the  solution  of  the  PDE  does  not 
take  its  maximum  in  the  Interior  of  the  domain. 

An  adaptive  mesh  technique  also  puts  the  additional  constraint  on  a  finite 
difference  scheme.  The  mesh  will  certainly  be  nonuniform  and  novine.  Thi= 
implies  that  a  difference  scheme  must  be  designed  for  a  moving  mesh.  To 
understand  how  to  generate  a  difference  scheme,  the  partial  differential  equation 
will  be  recast  in  a  moving  coordinate  system. 

Suppose  X  =  x(t)  is  a  point  moving  with  velocity  dx/dt  =  v  =  v  (x,t). 
Further  suppose  one  wants  to  follow  the  solution  of  (1.16)  as  the  point  x  moves. 
Then 

u(x,t)  =  u[x(t),t]   .  (1.36) 

u  is  now  a  function  only  on  time  so  that 
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du       3u   dx    .   3u  3u    ,    9u  .^    .-,. 

dI  =  3l^dr^^="8a^^3T    •  ^'-^'^ 


This     can    be    made  clearer    by    using    infinitesimals    as     in  Fig.    3.       By    Taylor's 
theorem: 

u[x(ti    +  dt),  ti    +  dt]    =  u[x(ti),    t, ]    +  dx  iZ    I        +  dt  —    I                         (1. 38) 

^  ^                                  ^          J-                 8x    It^  3t    It^ 


or 


du       3  u   dx  _  8  u 
dt        rx   dt   ~  3T 


Substituting  (1.38)  into  (1.16) 


2 
du  ,  ,       .  3u     3  u  .   _Q. 

dt         8  3x     3^2 


This  form  can  be  discretized  directly.   Before  doing  so  the  following  notation  is 
introduced.   Let 

j  =  1,  ...,  N 
Xj"  =  x.(nAt)  (I.^O) 

n  =  0,  1,  ... ,  M 


be  the  j*"   mesh  point  at  the  n^   time  level. 
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j  =  1,  ....  N 

n  =  0,  1,  ...  ,  M 

be  the  approximate  solution  at  the  j^"  point  and  time  level  n, 

j  =  1,  • • • ,  N 


X."  -  X  "-1 


v^ 


n 


J      J  (I.A2) 


J         At 

n  =  1 ,  . .  .  ,  M 


be  the  velocity  of  the  j*"   mesh  point  at  time  level  n, 


j  =  2,  . . . ,  N 
h.n  =  x.n  -  x._,n  (1.^3) 

"  =  C ,  .  ,  .... 


be  the  mesh  spacing  at  the  point  j  and  time  level  n, 


At  =  I   M  c  Z  (I-^^) 

M 


be  the  time  step  need  to  reach  a  time  T  in  M  steps. 

With  the  above  notation  the  following  consistent  finite  difference  schene  If 
proposed: 
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n+1  _  A  n  -  n+1  _  .    n+1 


At  J  h  "■'■^ 


X    n+li,  n+1  _  cv,  n+1  .  y,        n+l\A  n+1  .  x    n+lu   n+1 
*j+l   ^j      ^^j    ^  ^j+1    ^*j    ^  *j-l   ^j+1 


h, 


n-ri.    n-M/v,  nt-i  ■  v,   n-M-v/oi 


for  c  -  v,"+l  >  0 


r'-*:\<,_..n«,  "V^'-*r> 


At  J  V,   n+1 


.^^n+lh.n+1  _  (h.n+1  +  h.^^^^^)   *  ."+1  +  * -j_i"^^h  ^+^"+1 


< 


hj^-^^h.+i^+^Chj^+l  +  hj+i"+l)/2 


for  c  -  vx""^^  <  0 


^^"  =1   *n"  =  1   *j°  =  ^^'^j") 


(1.^5) 


(I.A6) 


(I.A7) 


This  scheme  preserves  the  maximum  principle  of  the  PDE.  To  see  that  the  maximum 
must  be  attained  at  the  boundary,  suppose  (^-i^  is  the  first  interior  maxima 
(1  <  3  <  N,  1  <  n  <  M)  in  the  sense  that 


n+1  >  ,|,  n+1  (j  <  3   ,   n  <  H)   .  (I.A8) 

3       ^ 
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Then  the  right  hand  side  of  (1.45)  or  (1.46)  is  positive  while  the  left  hand  side 
Is  negative;  a  contradiction. 

A  similar  argument  holds  for  a  minimum  value  being  attained  at  the  boundary. 
From  the  maximum  principle  it  is  seen  that  the  scheme  is  stable  in  the  maximum 
norm 


."U  =  max  l4)  ."!   .  (1.49) 


These  arguments  are  an  extension  of  ideas  described  in  Richtmyer  and  Morton  [IB]. 

With  the  stability  question  answered,  the  matter  of  convergence  can  be  dealt 

with. 

Theorem  3   Suppose  v^",  for  all  meshes,  is  bounded  everywhere  within  the  domain 

of  interest  [0,1]  x  [0,T].   Then,  (1.45)  to  (1.47)  is  a  convergent 

scheme  provided 


c  +  max  II  v"ll  -  ■ 

0  <  c^  <  ( ]  At  <  C2  (1.51) 

min  II  h"ll 
n 


max  h  J 

-^ <  C3   .  (1.52) 

min  hj" 


We  leave  the  proof  to  Appendix  B. 
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Asymptotlc  Theory 

Suppose  the  mesh  spacing  and  time  step  of  the  difference  scheme  in  the 
previous  section  are  so  small  that  all  the  quantities  in  (B.6)  can  be  replaced 
with  values  at  the  mesh  point  (x.",T^t)  without  significant  error.  Then  (R.6) 
becomes 


E"  =  ^  L?*  (x.",nit)  +  ^  (v.")2  !i*  (X  ",nAt)  -  v  "  1^   (x^",nAt) 


1  +  sgn(c  -  v.")  (c  -  vj")  hj"3  2. 

^        2  2     At  9^2  ^  J  '    ^ 


(1.52) 


l.,g,(,  _.n)  (c-vjn)hj,i«32,  ^^^^^^^ 


^^   3x2    ^ 


[h.,,^^^(x.n.  n.t)-h/L!*(x.n.  Mt)] 
3x  3x 


+  K 


3(hj+^"  +  hj")At 


Next,  assume  that  max  II  v"ll  <  2c  (where  c  is  the  advection  velocity),  in  (1.51) 

n 

C2  =  1,  and  adjacent  cells  are  nearly  equal  (h^"  =  h.^^").   Then 


-22- 

(1.53) 


+  h  2  (  I  ^  I  + 1^  I  !4  n 


3x2     9c   3^3 


In  this  form,  the  truncation  error  is  expressed  as  a  local  mesh  spacing  squared 
times  a  point  dependent  weighting  function.  This  is  just  the  form  of  the  error 
discussed  in  the  first  section  of  the  chapter.   By  setting 


w(x,t)  -J-[     I  i!*  I  +  Ac2  I  !^  I  +  c  I  1^  M 
ISc^     3t^  ^x^  ^'^^^ 


(I.5A) 


3x2     9^   3x3 


we  can  use  (I. A)  to  generate  a  mesh.   This,  of  course,  is  a  special  case  (t  =  1) 
of  (1.12). 

To  see  how  good  an  asymptotic  approximation  to  the  partial  differential 
equation  the  scheme  is,  it  is  necessary  to  examine  (1.8)  a  bit  closer.  Equation 
(1.8)  is  written  here  for  convenience 

wxf2  =  c^2  or  w^^2x^  =  ^^      ^  (1.8) 

2 

It  is  first  necessary  to  evaluate  c-^    .      This  will  tell  us  the  maximum  truncation 

error.  If  the  truncation  error  is  decreased  by  a  factor  k  then  the  error  will  be 
decreased  by  a  factor  k.  The  benchmark  to  which  the  above  is  compared  is  the 
truncation  error  for  the  uniform  mesh.  (Implicit  in  this  comparison  is  the  fact 
that  time  steps  used  in  the  adaptive  and  uniform  schemes  must  be  equal.) 
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max  w 


(N  -  1)2 

By  calculating  Cj^  in  (1.8)  the  truncation  error  for  the  nonuniform  mesh  can  be 
related  to  the  integral  of  the  square  root  of  the  normalized  weight  function. 

/  \        "^"   ^  )  dx  2 
ic^   max  w(C  )  -' 

ci2=„>ax!w|[  ^^-^^ ]    .  (1.55) 

Hence,  if  the  integral  of  the  square  root  of  the  normalized  weight  function  w  is 
small,  the  truncation  error  is  significantly  reduced.  By  examining  Fig.  2,  it 
follows   that  (1.55)  for  small  k  has  a  small  integral. 

It  may  now  be  verified  that  using  (1.4),  the  assumption  of  small  mesh 
variation  is  true.  Small  mesh  variation  would  be  indicate  by  a  small  value  of 
Xcr  .  Xrr  could  then  be  calculated  from  (1.8)  by  differentiation.  However  if  w 
is  not  continuously  dif f erentiable  we  cannot  use  the  mean  value  theorem.  One  can 
still  measure  the  mesh  variation  as  follows.   Take  two  adjacent  points  x.,x.. 


Then 


v(x^)^^(i)    =  c^^ 


The  time  derivatives  (fi  ^^   and  ({>  .^   can  be  expressed  via  the  PDE  in  terms  of  space 
derivatives  up  to  the  fourth  order. 
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(1.56) 


w(Xj)x^2(j)  =  ^^2 


Subtracting 


w(Xi)(x^(l)  -  x^(j))(x^(i)  +  x^(j))  +  (w(Xi)  -  w(xj))x52(j)  =  0   .    (1.57) 


Solving  for  the  mesh  spacing  difference 


2 
[w(xj^)  -  w(xj)]c^ 


(1.58) 


X^(l)-X^(j)     „(Xi)w(Xj)[x^(i)  +  X5(j)] 


[w(xj^)  -  w(Xj)]c 
[w(xi)w(Xj)]^/2f„l/2(^^)  ^  w^/2(x^)j 


Since  c  =  0(1/N-1)  the  mesh  variation  can  be  made  arbitrarily  small. 

From  experience  with  practical  computations  it  is  often  found  that  the 
weight  function  must  be  scaled  because  of  an  insignificant  number  of  mesh  points. 
One  can  then  ask  how  much  (1.55)  is  degraded  by  a  scaling  of  the  weight  function. 
Let 


-*  =  -(-  -  -min>/(-max  "  "min^  ^  '  (^-^9) 


where 


"max  =  '^'^  "('^^    •    "min  =  ""i"  "('^> 
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and 


1  <  a  < 


w     ~  w  • 

max    min 


"min 


Solving  (1.8)  for  c^ 


N   *l/2  , 
*  2     2   f     r        w  ^' ^  dx  >, 

^1 


(1.60) 


using  (1.59) 


2    max   *min   ?   r   max   *niin       > 


a 


\' 


w    —  w  • 
max    min   2 
<  Ci 


Letting 


(1.61) 


a  =  9 


w     ~  w  « 

max    min 


min 


w. 


min 


max    min 


<  e  <  i) 


(1.62) 


(1.61)  becomes 


w_ 


1/2 


wx.   =  [  , dxj 

x,  (N  -  \y 


(1.63) 
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For  e  =  1  the  formula  reduces  to  the  original  formula  with  weight  function  w.  The 

2 
following  estimate  will  be  made  to  give  a  crude  idea  how  wxr   behaves  away  from 

e  =  1. 


1/2  1/2 


-  wl/2 


(I.6M 

w  .   f  i-lll 
=  „l/2  +  min  ^    fi    J 


[-^-mlnt  ^)]'^'^wl/2 


<  „i/2  ^  '"^"^   e   ^ 


1/2 


2w 


^     2   w   '^    6    ^-l 


<  wi/Mi-i[ii^)] 


The  factor,  for  6  near  1 


[l+i[-i-f^)]  (1-65) 


gives  the  increase  in  truncation  error  as  a  result  of  rescaling.  From  nunerical 
experience  the  estimate  is  overly  pessimistic  as  shall  be  seen  in  the  following 
sections. 
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Comblned  Algorithm 

The  asymptotic  theory  of  the  last  section  gives  only  a  rough  guide  to  the 
performance  of  a  practical  algorithm.  In  fact,  the  regime  of  Interest  Is  the  one 
where  as  few  points  as  possible  are  used.  However,  the  previous  theory  can  still 
be  helpful  in  constructing  a  finite  difference  scheme.  The  rest  of  this  section 
is  devoted  to  the  description  of  a  practical  scheme. 

To  describe  a  practical  scheme  It  is  best  to  start  with  a  general 
description  of  the  algorithm.  Suppose  at  time  level  (n)  one  has  data  on  a  mesh. 
To  solve  (1.45,47),  a  mesh  at  time  level  (n+1)  is  required  along  with  the 
corresponding  mesh  velocities.  To  generate  a  mesh  at  time  level  (n+1)  a  weight 
function  w  is  generated  at  the  time  level  n  from  the  given  data.  Then  (1.15)  is 
used  to  calculate  a  new  set  of  points 


Xj    j  =  1,  2,  ...,  N   .  (1.66) 


Next,  the  quantity 


Ix   -  Xj"l 

%ax  =  ^^^    ^  ^,  (1-67) 


is  calculated.   Finally  if  v    is  less  than  the  given  bound  v.  (a  bound  on  the 
mesh  velocity  needed  for  convergence) 


Xj""*"!  =  Xj  (1.68) 


otherwise 


^"^1  =  xj"  +  (Xj  -  Xj")  vi,/v^^^   .  (1.69) 
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Using  (1,45-47),  the  data  may  be  advanced  to  time  level  (n+1)  and  the  above 
cycle  may  be  repeated.   This  algorithm  is  summarized  in  figure  A. 

We  now  fill  in  some  of  the  details  of  the  flow  chart.  To  calculate  a  new 
mesh  equation  (1.15)  is  used.  Equation  (1.15)  is  a  two  point  boundary  value 
problem.  As  a  result  of  the  type  of  problem  (1.15)  is,  one  has  a  treTnendous 
variety  of  methods  to  choose  from  (see  [6B]).  The  method  chosen  here  is  a  finite 
difference  method.        '      '   ;    ■ 


••  '    _i    [w(Xj+i)  -  w(Xj_^)](Xj+i  ~  ^j-l) 

"'^i+l  -  2xj  +  Xj_i  =  —  t  ^ 

■^  J    J      8    j8(i  _  t)  +  (Xj+;^  -  *Xj_i)^w(xj)t] 


(1.70) 


j  =  2,  ...,  N-1   . 

The  system  of  equations  can  be  solved  iteratively 

(m)       (m)     (m)     (n)  4 
(m+1)      (m+1)     (m+1)    _i      M'^  j+0    "  w(x  ^.^ )  ]  (x  .^^  -  x^.^) 
Xj+i    -  2x.  +   Xj_i    =     t  

(m)     (m)  3    (m) 
[8(1  -  t)  +  (x^^^    -   x._^)   w(Xj   )t] 


■  A  !   -1-3 


where  the  superscript  (m)  is  the  m*"   iterate.   A  relaxation  parameter  X  is  also 
used: 
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(m)       (m)    (m)    (m)  4 
(m+1)           (nH-1)     (m+1)    _i     l^^C^^+O    "  "(x  j_i ) )  (Xj+j  -  x._^) 
Xj+i   -  (2  +  X)xj     +  x^_i  '  =  J:  t 1 :' (1.72) 

(m)    (m)  3    (m) 
[8(1  -  t)  +  (xj+i  -  x._^)   w(Xj   )t] 


AXj 


Finally  the  solution  Is  rescaled  so  that  monotonicity  is  preserved.   Vhen  the 
following  inequality  is  satisfied 

,  (m+1)     (m),  (m)     (m)    (m)    (m) 

Ixj     -  Xj   I  <  Y^iin  [*Xj    -  Xj_j,  Xj^^  -  x^   ],   0  <  Y  <  1 


then 


(m+1)     (m+1) 
xj      =xj    .  ^.  (1.73) 


However  when  the  inequality  is  not  satisfied 

(m+1)     (m)  (m)   Jm)     Jm)        Jm) 

Ixj     -  Xj   1  >  Ymin  [Xj    -  Xj_p  Xj^^  -  x^   ]  (I.7A) 


then 


(m+l)      (m) 
(m+1)     (m+1)    X.     -  X.  (m)     (m)    (m)     (m) 


This  iterative  procedure  is  chosen  because  it  uses  little  computation  time.   The 
LU  factors  for  the  left  hand  side  need  be  calculated  only  once.   Since  the  scheme 
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is  used  at  every  time  step  of  a  calculation  to  generate  a  new  mesh,  the  old  raesh 
Is  used  as  a  first  guess  for  the  iteration.   If  the  weight  function  does  not 
change  drastically  from  the  previous  one,  then  the  old  mesh  is  an  accurate  first 
guess  and  only  a  few  iterations  are  required  to  generate  a  new  raesh. 

The  last  part  of  the  algorithm  to  describe  is  the  weight  function.  *Caively 
one  could  use  the  PDE  (1.16)  to  convert  all  the  terms  with  time  derivatives  in 
(1.54)  to  space  derivatives  and  then  use  the  new  form  as  a  weight  function. 
However,  since  few  points  are  used,  the  higher  order  differences  will  be  both 
costly  to  compute  and  noisy.  Instead,  just  the  lowest  order  truncation  error 
terras  will  be  used. 


w(x,t)  =  [    *  )   .  (1.75) 

3x^ 


Independent  of  the  truncation  error  one  can  give  justification  to  the  above 
weight  function  by  the  following  argument.  The  local  error  made  by  a  piecevoise 
linear  approximation  to  a  twice  dif f erentiable  function  is  proportional  to  the 
second  derivative  at  the  function  times  the  mesh  spacing  squared.  By  placing 
more  mesh  points  in  an  area  of  high  curvature,  a  better  approximation  can  be 
made.  Since  the  approximate  solution  to  the  PDE  is  piecewise  linear  it  would 
best  represent  data  by  bunching  points  in  high  curvature  regions. 

Another  candidate  for  use  as  a  weighting  function  might  be  the  gradient 
length  function  used  in  [9A].   '       '  ■•'   ., 


*    2 
w(x,t)  =  [  -^  )  (1.76) 


An  argument  for  its  use  is  as  follows.   Consider 
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e'^^  (1.77) 


then 


w(x,t)  =  6^  (I.7P) 


when  6  is  large,  the  function  ^    would  require  many  points  to  accurately  resolve 
it.   When  6  is  small  fewer  points  are  required. 

One  additional  step  that  is  needed  to  complete  the  discussion  on  weight 
functions  is  smoothing  and  scaling.  If  a  small  nuraher  of  points  is  used  with  a 
weight  function  with  large  variation,  the  large  variation  in  adjacent  mesh  widths 
will  cause  a  loss  of  accuracy  in  the  differencing.  Computational  experience  has 
also  shown  a  poor  rate  of  convergence  when  the  weight  function  becomes  more 
jagged  or  becomes  larger  in  variation.  A  practical  smoothing  method  is  the 
following  iterative  scheme. 

„(m)(5^^n)  =  „(m-l)(-^^n)  +  g[  (^(m-l)(^___^^n)  (1.79) 


-  Zw^'^'^^Xj")  +  w^'^-l^Cxj.^")] 


where  3  ^  111  for  reasons  of  stability.  (This  is  the  explicit  approximation  for 
the  heat  equation  with  mesh  and  time  spacing  1  [IB].)  In  addition  to  smoothing 
the  weight  function,  the  "diffusion"  of  the  weight  function  hunches  the  mesh  in  a 
region  that  borders  where  the  weight  function  is  large.  This  diffusion  seems  to 
compensate  for  not  using  all  the  truncation  errors  terms  (compare  figures  in  the 
second  section  to  see  how  the  higher  derivatives  extend  further  from  the  boundary 
layer).   Scaling  was  Introduced  in  the  section  on  asymptotic  theory. 
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Numertcal  Experiments 

In  this  section,  two  numerical  experiments  are  described.  The  first 
experiment  will  deal  with  the  size  of  the  time  step  while  the  second  experiment 
will  deal  with  the  number  of  mesh  points. 

First  Experiment 

Since  the  finite  difference  scheme  described  in  the  previous  section  is 
unconditionally  stable,  any  time  step  may  be  used.  However,  as  will  be  seen  from 
experiment,  accuracy  conditions  give  a  bound  on  how  large  a  time  step  should  be 
used. 

Instead  of  using  "At"  as  a  description  of  the  time  step,  it  will  be  nore 
convenient  to  use  a  dimensionless  quantity  that  has  the  forra  of  a  Courant  number 


(c  +  Iv,  I)  At 

- .  (1.80) 

Ax 


Here  c  is  the  advection  velocity,  v,  the  bound  on  the  mesh  velocity,  and  Ax  the 
minimum  mesh  spacing.  v^  will  be  zero  and  Ax  inversely  proportional  to  the 
number  of  points  for  the  uniform  mesh. 

From  the  above  formula,  it  is  easy  to  see  that  a  nonuniform  adaptive  schenc 
will  always  have  a  smaller  timestep  than  a  uniform  fixed  mesh  for  a  fixed  a.      The 
first  experiment  will  use  the  following  values, 
c  =  1  (c  is  the  advection  velocity) 

<    =  0.01  (k  is  the  diffusion  constant) 

NPTS  =  51  (NPTS  is  the  number  of  mesh  points) 

X  =  1  (X  is  the  relaxation  parameter  in  1.72) 
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t  =  0.9,0  (t  is  the  interpolation  parameter  in  1.14  with 

0  for  uniform  mesh  and 

0.9  for  adaptive  mesh) 
Y  =  0.8  (Y  is  the  mesh  excursion  liraiter  in  1.73) 

''  ~  "max'^^min  ~  ^'^'^   ^^  ^^     ^^^  ratio  of  the  maximum  to  minimum  values  of  the 

weight  function) 
nsmth  =  20  (nsmth  is  the  number  of  smoothing  iterates  done  in  1.79) 

niter  =  3,100       (niter  is  the  number  of  solution  iterates  to  1.73) 

100  for  the  initialization  of  the  mesh  and  3  for  each 

succeeding  time  step 
6  =  0.4  (6  is  the  smoothing  parameter  in  1.79) 

TFIN  =  0.1  (TFIN  is  the  ending  time  for  the  calculation) 

Iv^l  =  4,0  (l^b'  ^^  ^^^  upperbound  allowed  for  the  mesh  velocity  in 

1.69  with 

0  for  the  uniform  mesh  and 

4  for  the  nonuniform  mesh) 
A  =  0.01  (A  is  the  initial  boundary  layer  thickness  in  1.34) 

The  weight  function  chosen  is  the  curvature  8  <^ /3  x  and  is  calculated  using  the 
differencing  of  the  right  hand  side  of  equation  (1.45).  The  curvature  is  then 
rescaled  so  that  r  =  w^^^  =  100  and  w^j^^^^  =1.  In  the  experiment  a  will  vary 
between  10  and  10"  in  logarithmic  increments  of  10^'  .  Table  1  and  Fig.  5 
summarize  the  results  of  the  experiment.  It  should  be  emphasized  that  for  a 
given  Courant  number  the  time  steps  for  a  uniform  and  nonuniform  mesh  will  differ 
because  of  different  mesh  spacings  and  mesh  velocities. 

There  are  several  comments  about  the  results  that  can  be  made.  First,  it 
should  be  noted  that  the  asymptotic  error  differs  by  slightly  more  than  a  factor 
of  ten. 
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The  Eq.  (1.8)  gives  a  factor  of  ten  difference  in  mesh  spacing  which  consequently 
reduces  the  truncation  error  by  this  factor  in  the  region  where  the  curvature  is 
large.  At  the  same  time  the  mesh  spacing  isn't  so  large  in  the  region  where  the 
curvature  is  small  so  as  to  enhance  the  truncation  errors  there  significantly. 
Thus  overall  the  maximum  truncation  error  is  reduced  by  a  factor  of  ten  and  hence 
the  absolute  error. 

A  second  comment  concerns  the  Courant  number's  effect  on  the  mesh  generator. 
As  the  Courant  number  goes  to  zero,  the  mesh  generator  has  more  time  to  adapt  to 
the  problem.  This  is  indicated  by  the  stabilization  of  the  mesh  ratio,  in 
Table  1 ,  as  a  goes  to  0.01. 

Figures  6  thru  10  show  particular  solutions  and  weight  functions  taken  from 
runs  used  in  generating  Table  1.   Figure  6  is  an  initial  configuration. 

The  grid  mapping  can  be  used  to  see  where  the  mesh  points  are  by  the 
following  procedure.  First,  superimpose  a  uniform  mesh  on  the  Y-AXIS.  Second 
draw  lines  from  the  mesh  points  on  the  Y-AXIS,  parallel  to  the  X-AXIS,  to  the 
curve.  Finally,  from  the  points  of  intersection  of  the  lines  and  the  curve,  draw 
lines  parallel  to  the  Y-AXIS  to  the  X-AXIS.  The  points  of  intersection  on  the 
X-AXIS  are  the  locations  of  the  mesh  points.  The  most  important  feature  to 
notice  is  how  well  resolved  the  "corners"  of  the  initial  data  are.  The  "corners" 
are  of  course  the  places  where  the  curvature  is  highest. 

Figure  7  shows  simultaneously  the  computed  solution,  absolute  error,  crid 
mapping,  and  analytic  solution  at  time  t  =  1.  It  should  be  noticed  that  the 
error  is  concentrated  in  the  regions  of  large  curvature. 

Figure  8  shows  the  scaled  and  smoothed  weight  function.  Notice  the  width  of 
the  weight  function  is  large.  This  will  cause  adjacent  mesh  spacings  around  the 
high  curvature  region  to  be  small  as  well.  This  point  is  crucial  for  two 
reasons.   The  first  reason  is  that  the  "disturbance"  or  "front"  is  moving  with 
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the  advection  velocity  c.  Since  the  mesh  generator  uses  the  n  time  level  for 
generating  a  mesh,  spreading  the  weight  function  ensures  that  the  disturbance 
will  not  move  out  of  a  region  of  high  resolution  in  a  time  step.  Thus  the 
smoothed  weight  function  has  a  safety  margin  built  into  it.  This  safety  margin 
is  also  important  for  another  reason.  By  examining  the  figures  in  the  second 
section  it  is  fairly  obvious  that  higher  derivatives,  which  also  contribute  to 
the  truncation  error,  have  maxima  and  minima  in  different  spots  than  does  the 
curvature.  From  experiment,  smoothing  or  smearing  the  weight  function  minimizes 
higher  order  truncation  errors  since  all  truncation  errors  are  large  only  in  a 
neighborhood  of  the  wave  front. 

The  9th  and  10th  figures  are  the  corresponding  solutions  on  a  uniform  mesh. 
The  absolute  error  is  similar  to  the  nonuniform  case  except  a  factor  of  ten 
larger.  Also  apparent  is  less  resolution  of  the  initial  data  as  well  as  a  more 
diffuse  solution  profile. 

Second  Experiment 

The  following  experiment  will  hold  the  Courant  number  fixed  and  vary  the 
number  of  mesh  points  used.  Here,  another  dimensionless  constant  will  be  useful. 
Define 


&       < 


This  will  be  called  the  grid  Reynolds  number  because  of  similar  structure  to  a 
fluid  Reynolds  number.  Ax  is  defined  as 


Ax  =  1 .  (1.82) 

npts  -  1 
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where  npts  Is  the  number  of  points  used  in  the  computation  mesh.   In  this  case  Ax 
has  meaning  for  both  uniform  and  nonuniform  meshes. 

The  experiment  will  fix  the  Courant  number  a  at  0.1,  <  and  L  will  be 
decreased  to  0.005  and  the  smoothing  Iterates  will  depend  quadratically  on  the 
number  of  points  used 

2 

nsmth  =  (  nZll]         .  (1.83) 

^  11.4  '' 

This  insures  approximately  the  same  weight  function  independent  of  the  number  of 
points.  Notice  that  when  npts  =  51  nsmth  =  20.  For  npts  large  nsmth  will  be 
impractical  for  real  problems.  This  can  be  circumvented  by  using  an  implicit 
algorithm  for  the  smoothing.  There  would  be  no  real  increase  in  cost  of 
smoothing  since  the  LU  factors  are  already  stored  for  use  with  the  mesh  generator 
algorithm.   Table  2  and  Fig.  11  give  the  results  of  the  experiment. 

There  are  several  observations  to  be  made  about  the  results  of  the  numerical 
experiment.  Perhaps  the  most  striking  observation  is  the  constancy  of  the  error 
over  a  wide  range  of  values  of  R  for  the  adaptive  mesh.  A  similar  explanation 
to  the  one  given  in  the  first  experiment  tells  us  why  this  Is  so.  The  adaotive 
mesh  is  minimizing  the  large  truncation  errors  in  regions  of  strong  variation 
while  not  increasing  the  truncation  errors  significantly  elsewhere. 

Another  observation  is  the  differences  in  computation  time.  The  differences 
are  brought  about  by  a  Courant  number  imposed  on  the  problem.  In  all  cases  the 
adaptive  scheme  does  better.  However,  what  has  to  be  weighted  is  the  increased 
time  used. 

For  R  <  1  the  errors  come  close  together  in  an  absolute  sense.  This  is 
because  the  uniform  mesh  starts  to  resolve  the  problem.  But  the  real  utility  of 
the  adaptive  scheme  is  that  It  extends  the  accuracy  of  a  difference  method  for 
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R  >  1  making  it  Indispensible  in  problems  where  the  user  has  a  limited  number  of 
points  or  in  very  singular  problems.      ■..••■   ;  .   •  ,■  ;  ;>•         .  • 


(TABLE  2) 


npcs 

'& 

absolute 

error 

1 

nsmth 

mesh  ratio 

time  steps 

uniform 

nonuniform 

uniform 

nonunl form 

801 

1/4 

.1524  X  lO'J 

.8246  X  10"^ 

4937 

13.7 

801 

8498 

401 

1/2 

.2881  »  10"' 

.1033  X  10"| 

1237 

15.8 

401 

5542 

268 

3/4 

.4092  K  10"' 

.1308  X  10"' 

553 

14.4 

268 

3778 

201 

1 

.5211  X  10"' 

.1533  X  10"' 

311 

13.4 

201 

2899 

161 

1  1/4 

.6218  X  10"' 

.1714  X  10"' 

199 

12.8 

161 

2387 

134 

1  1/2 

.7194  X  10"' 

.1855  X  10"' 

138 

12.3 

134 

2051 

115 

1  3/4 

.8032  X  10"' 

.1956  X  10"' 

102 

11.9 

115 

1819 

101 

2 

.8801  X  10"' 

.2021  X  10"' 

78 

11.6 

101 

1655 

90 

2  1/4 

.9618  X  10"' 

.2026  X  10"' 

62 

11.3 

90 

1525 

81 

2  1/2 

.1048 

.1980  X  10"' 

50 

11.1 

81 

1419 

74 

2  3/4 

.1103 

.1879  X  10"' 

42 

11.0 

74 

1331 

68 

3 

.1162 

.1723  X  10"' 

36 

10.8 

68 

1249 

63 

3  1/4 

.1249 

.1542  X  10"' 

31 

10.8 

63 

1180 

58 

3  1/2 

.1297 

.1322  X  10"' 

26 

10.7 

58 

1110 

54 

3  3/4 

.1351 

.1105  X  10  ' 

22 

10.7 

54 

1055 

'' 

4 

.1434 

.9884  «  10"^ 

20 

10.7 

51 

1003 
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II.   2-D  Grid  Generation 

All  the  Ideas  about  mesh  generation  of  the  previous  chapter  can  be 
generalized  to  higher  dimensions.  In  fact,  as  will  be  seen,  higher  dimensions 
have  possibilities  for  variational  principles  that  are  not  seen  in  a  single 
dimension.   To  start  with  the  parameter  and  physical  spaces  will  be  introduced. 

The  parameter  space,  which  is  really  a  set  of  array  indices,  will  now  have 
two  variables.  Instead  of  using  greek  letters,  the  parameter  variables  will  be 
u,v  (see  Fig.  12).  The  parameter  space  will  be  the  rectangle  where  IMAX  and  JMAX 
following  fortran  conventions  are  integers. 

The  physical  space  can  be  any  multiply  connected  region  in  two  dimensions 
that  can  be  reduced  by  cutting  to  a  simply  connected  region. 

The  numbers  in  figure  13  indicate  orientation  and  boundries.  As  shown,  the 
variable  names  for  the  physical  space  are  x,y. 

With  the  spaces  defined,  it  is  now  possible  to  turn  to  the  variational 
formulations.   The  first  functional  to  be  described  is  the  smoothness  functional 


n 


.....  (ii.i) 

=  f   u'^+u^  +  v^+v^dX 
J    X     y     X     y 

where  dX  is  the  physical  space  volume  element  and  fJ  is  the  domain  of  integration 
in  the  physical  space.   The  Euler  equations  in  physical  coordinates  are, 


^5^*l^5^<"^<'*V*"x^*v/,  -   2(„,,.„^p.O         (11.2) 
X       y 
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i_  J_  +  i_  A-  (u^2  +   2  +   2  +   2)  ^   2(v^^  +  v   )=0   .      (II. 3) 
X        y 


Using  the  maximum  principle,  it  is  easy  to  see  that  all  the  coordinate  lines 
(u  =  const,  or  V  =  const.)  are  within  the  physical  region  and  that  the  solutions 
to  the  Laplace  equations  minimize  the  variation.  This  is  of  course,  the  2-D 
analogue  of  the  1-D  smoothness  functional  in  the  previous  chapter. 

The  second  functional  to  introduce  is  the  weighting  functional.   By  analogy 
to  the  first  chapter 


I„(u.v)  =  /   w(x,y)  J  dX  (II. A) 


where 


J  =  det  I  ''^  ''M   .  ■  (II-5) 

^u  yv 

Notice  here  that  instead  of  controlling  mesh  spacing  directly  the  variational 
principle  controls  mesh  spacing  by  controlling  the  size  of  area  elements. 

There  are  advantages  and  disadvantages  to  this  principle.  A  disadvantage  to 
the  principle  is  the  fact  that  the  Euler  equations  for  the  variational  principle 
have  many  solutions  for  a  given  weighting  function.  This  has  been  shown  in  an 
earlier  paper  [9A].   As  a  result  the  computation  is  certainly  not  well  posed. 

A.  M.  Winslow,  [29A],  has  suggested  a  different  variational  principle  to 
circumvent  the  problem  of  nonuniqueness,  written 

f,i  .      .  „  ■■     ..■■■■■',■■' 

Iw'  =  /  «(^.y)  IC^^)"^   +  (^n)2]  dX   .  (II. 6) 
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The  corresponding  Euler  equations   are   linear,   uncoupled,   and  have  unique 
solutions.   ■".■•.■.'     ,   . 

However,  when  the  problem  of  nonunlqueness  can  be  overcome  the  variational 
problem  I^  is  more  useful.   This  stems  from  the  fact  that  the  property 

wJ^  =  const.  (II.7) 

can  be  recovered  from  the  Euler  equations.   This  formula  is  quite  useful  in 
scaling  the  weight  function  property  as  discussed  in  the  one  dimensional  section. 
To  overcome  the  nonunlqueness  problem  the  functionals  I   and  I   are  added 
together 


I  =^sls  +^w^w   •,   .....  ...,.  (II-8) 


It  is  conjectured  that  the  smoothing  functional  combined  with  the  weighting 

functional  yields  Euler  equations  which  have  a  unique  solution.   Further,  it  is 

2 
conjectured  that  as  X^  gets  large,  the  property  wJ  =  const.   is  approached. 

Numerical  calculations  to  support  these  conjectures  will  be  presented. 

Another  variational  principle  that  will  prove  useful  is  a  functional  that 

measures  orthogonality  of  mesh  lines 


Iq  =  /  [07u)  •  (Vv)]2  j3  dX   .  (II. 9) 

The  J   term  weights  the  measure  in  favor  of  orthogonalizing  regions  of  larger 
cell  areas.   The  more  obvious  functional  would  be 

^o  =  /  t^")  •  C^v)]^  ^^      '  (11.10) 
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Both  functlonals  yield  Euler  equations  that  don't  have  unique  solutions  as  well. 
Again,  to  circumvent  this  difficulty  I^  or  I   will  be  used  in  conjunction  with 
the  smoothing  function.   So  now  the  variational  principle  looks  like 

I  =  Xg  /  [(Vu)2  +  (Vv)2]  dX  +  X^  /  wJdX  +  Xq  /  (Vu  •  Vv)^  J^dX   .     (II. U) 

So  where  one  combines  a  set  of  functionals  for  the  one-dimensional  case  as  a 
matter  of  convenience,  in  higher  dimensions  one  combines  the  set  of  functionals 

as  a  matter  of  necessity.  •-,         -;    ,  "  ■ 

Normalization  and  change  of  variables  ' 

Before  doing  any  computations  it  would  be  advantageous  to  prepare  the 
variational  principle  in  several  ways.  The  first  change  would  be  to  normalize 
the  Integrals,  I^^  and  I^,  so  they  will  be  of  the  order  of  the  smoothing  integral. 
Normalization  factors  are  introduced. 


I  =  _f_  /  [(Vu)^  +  (Vv)2]  dX  +  -^  /  (Vu  •  Vv)^J^  dX  +  _!:^  /  wJ  dX   .  (11.12) 
^s'  ^o'  -^w 


A  fairly  effective  and  simple  way  to  determine  X^',  X^',  and  X^'  is  to  use 
dimensional  analysis.  First  X^'  is  set  to  one.  The  integrand  of  the  smoothing 
integral  has  dimensions  of  pararmeter  length  £  over  physical  length  i  squared. 
That  is 


~   2  •■   "■•■■■ 

I  =  [  -  )    .  (11.13) 
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Similarly  the  integrands  of  the  orthogonality  and  weighting  functionals  have 
dimensions        ■'■   ~ 

1^=  [  _]  (II. lA) 

I 


2 

Iv  =  w  [  i  )  (11.15) 


where  w  is  the  unit  of  measure  of  the  weight  function.   In  order  that  I   and  I 

°  o  V 

have    the   same   units   as   I  -;       -  ■•  ■    -.--        .  .    ! '  ■   <^     ,     ■•  -  ■-        •■  ■ 

A 

Aq'  =    (    ^)  (11.16) 


A 
^v'  =  ^  (    -  )         •  (11.17) 


Actual  values  used  in  calculations  are  the  following.  £  is  the  maximum  length  in 
the  physical  domain,  £  is  the  number  of  points  in  one  direction,  and  w  the  volune 
average  of  w. 


w  =  _  /  w  dX 
V  '' 


(11.18) 


where 


V  =  /  dX   .  (11.19) 
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Next  the  Independent  and  dependent  variables  are  interchanged  so  that  centered 
differences  can  be  used  in  approximating  the  various  derivatives  in  the  Euler 
equations.   The  variational  principle  becomes 


_  ,  r    (Vx)2  +  (Vy)2  ...  .  ^o  .  .      .      .2 


^o " 


(11.20) 


+  ^  j    w[x(u,v),y(u,v)]j2  dU 


where  dU  is  the  volume  element  in  parameter  space. 

Computations 

With  the  variational  principle  derived  the  Euler  equations  may  be  calculated 
from 


3x   3u  3x.,   3v  3x„  ^  s  i.        j       ) 


(11.21) 


o  w 


and 


_3_   _3_  3     i_  J_  [x   f  (Vx)^  +  (7y)^ 
3y   3u  3Xy   3v  3x^  ^  s  ».       j 
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CII.22) 


o  0    ^w    2^ 


Details  of  the  numerical  solution  of  the  above  problem  can  be  found  in  [9A]  and 
will  not  be  repeated  here. 

What  will  be  looked  at  are  two  particular  calculations.  The  first 
calculation  gives  evidence  that  the  combined  smoothing  and  weighting  functionals 
give  rise  to  equations  that  for  ^„  >>  ^s  bas  the  property  that 

wJ^  =  const.  (11.23) 

The  second  calculation  will  show  that  for  certain  classes  of  geometries  the 
variational  principle  will  not  generate  a  satisfactory  computation  mesh.  For 
these  geometries,  alternative  variational  principles  will  be  introduced. 

First  Example 

The  first  calculation  will  map  the  parameter  space  rectangle  into  a  unit 
square  in  the  physical  space  (shown  in  figure  14).  IMAX  and  JMAX  are  the  number 
(consistent  with  Fortran  notation)  of  points  used  in  each  direction.  The  total 
number  of  points  is  then  the  product  IMAXxJMAX.  The  numbers  1,  2,  3,  and  4 
indicate  the  relative  orientations  and  boundary  conditions. 

In  this  numerical  experiment  X^  Is  set  to  zero,  Xg  is  set  to  1  and  X^,  will 
vary.   The  weight  function  has  the  form 

w(x,y)  =  a  e  '^  '^  +  1  (11.24) 
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where 


r^  =  (x  -  1/2)2  +  (y  -  1/2)2  (11.25) 


and  a  and  s  are  positive  constants.   Table  3  shows  the  mean  over  the  variance  of 

2 

the  values  wJ   at  all  the  interior  mesh  points  as  X,,  varies.   Note  that  the 

w 

variance  diminishes  as  X^  increases  indicating  the  approach  to  the  solution 

9 

wJ"  =  const.   The  graph  in  figure  15  summarizes  the  results  of  the  table  3. 

Finally  the  following  contour  plots  in  figures  16  and  17  show  the  function  value 

wJ-  for  X,,  =  0  and  X,,  =  10.  ' 
w  w 


(TABLE  3) 

a  =  256 

s  =  1 

IMAX  =  6A 

JMAX  =  64 

^w 

variance/mean 

.001 

8.6A  X 

10-2 

.00316 

6.92  X 

10-2 

.01 

A. 53  X 

10-2 

.0316 

2.35  X 

10-2 

.1 

1.01  X 

10-2 

.316 

A. 27  X 

10-3 

1. 

2.15  X 

10-3 

3.16 

1.71  X 

10-3 

10. 

1.76  X 

10-3 
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Notice  all   the  gradients  are  at   the  boundary  where  the  points  are 
constrained. 

Second  Example 

The  second  example  will  have  the  parameter  space  mapped  into  a  circle  with  a 
slit  as  illustrated  in  Fig.  18.  The  arrows  and  numbers  indicate  orientation  and 
boundary  placement.  Notice  the  line  31  is  collapsed  to  a  point  in  the  physical 
plane. 

One  desires  that  the  mesh  look  like  a  discreet  version  of  polar  coordinates. 
However,  if  only  the  smoothing  functional  is  used,  lines  are  pulled  in  toward  the 
center  of  the  circle  and  eventually  converge  as  shown  in  Fig.  19. 

With  periodic  conditions,  as  shown  in  Fig.  20  at  an  intermediate  iteration. 
the  circles  are  collapsing  to  a  point  at  the  center  as  noted  earlier  by  Thompson 
[2A]. 

The  reason  for  the  problem  is  that  a  singularity  is  present  at  the  center  of 
the  circle.  That  is  u  is  a  logarithmic  function  of  the  distance  from  the  center. 
This  can  be  seen  by  solving  Laplace's  equation  in  cylindrical  coordinates. 

Remedy 

As  noted  by  Thompson,  a  radially  weighted  smoothing  functional  can  be 
derived  to  remedy  the  situation. 


l'=—-± _± +  —  -J1 ^1_  dU  (11.26) 


where  R^  is  a  scale  length  (to  preserve  units)  and 
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r  =  (x^  +  y2)l/2   .  (11.27) 


In  physical  space  the  functional  becomes 


Ig'  =  /  -L  (Vv)2  +  _^  (Vu)2  dX   .  (II. 2S) 


Using  the  cylindrical  volume  element 


I^'  =  /  —  r^(Vv)2  drde  +  /  R„  (Vu)^  drdB    ,  (11.29) 

•^o 


When  V  =  v(e  )  and  u  =  u(r)  the  functional  yields  the  uniform  spacing  desired. 
The  following  figures  demonstrate  the  remedy  for  a  circle.  Figure  21  shows  a 
polar  mesh  whose  interior  mesh  points  have  been  randomly  displaced.  The  boundary 
slit  is  made  periodic.  Figure  22  shows  the  mesh  after  20  iterations  with  the 
smoothing  algorithm  Ig'. 

The  weighting  functional  for  cylindrical  coordinates 

In  the  previous  section,  the  mesh  is  made  to  behave  in  a  one  dimensional 
manner  in  the  radial  direction  or  theta  direction.  By  modifying  the  weighting 
functional,  similar  behavior  can  be  obtained.   Let 


^w'  -  /  V^4  '^'^      '  (11.30) 

r 


2 
Where  R^   is  introduced  to  preserve  units.   Using  the  cylindrical  volume  element 

and  the  associated  coordinates 
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„'  =  /  Rq^w  i   drde   .  (11.31) 


I, 


When  V  =  v(9 )  and  u  =  u(r)  the  functional  reduces  to  the  1-D  case  since  the 
Jacobian  is  of  the  order  of  r  as  r  goes  to  zero- 
Figure  23  shows  a  solution  of  a  circular  mesh  problem  with  a  radially 
dependent  weight  function.  The  radial  distances  of  the  mesh  lines  are  then  shown 
in  the  following  table  along  with  values  that  were  generated  by  the  1-P  solver 
used  in  the  1-D  chapter. 
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( TABLE  4) 

Radial  Coordinate  Line  Locations 


2 
weight  function:   w(r)  =  255  e"^"*"-^^  ''^  +  1 


1-D         2-D 


0. 

0. 

0.10856 

0.10855 

0.17603 

0.17602 

0.22672 

0.22670 

0.26844 

0.26843 

0.30469 

0.30468 

0.33732 

0.33731 

0.36742 

0.36741 

0.39573 

0.39573 

0.42279 

0,42278 

0.44897 

0.44897 

0.47462 

0.47462 

0.50000 

0.50000 

0.52538 

0.52538 

0.55103 

0.55103 

0.57721 

0.57722 

0.60426 

0.60427 

0.63258 

0.63259 

0.66268 

0.66269 

0.69531 

0.69532 

0.73156 

0.73157 

0.77329 

0.77330 

0.82397 

0.82398 

0.89144 

0.89115 

1.0 

1.0 
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III.   Three  Dimensional  Grid  Generation 

Generation  of  grids  in  three  space  dimensions  is  a  direct  extension  of  the 
two  dimensional  case.  The  physical  coordinates  for  three  dimensions  are  x,  y, 
and  z  and  the  parameter  space  coordinates  (illustrated  in  figure  2A)  for  three 
dimensions  are  u,  v,  and  w.  However,  it  will  be  convenient  to  express  the 
physical  and  parameter  space  coordinates  using  superscripts  respectively  as  x 
and  u  . 

Derivatives  of  the  physical  coordinates  with  respect  to  the  parameter 
coordinates  and  vice  versa  will  be  denoted  by 


i   3x^ 
x^   =  r 


i,j  =  1,  2,  3  (III.l) 


and 


u.i  =i^  ,    i.j  =  1,  2,  3  (III. 2) 

^    3xJ 


So   that,    for  example, 


X32=?2  (III. 3) 

^         3w 


and 


3   =  i2i      .  (III. A) 

9x 
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Hlgher  derivatives  will  have  more  subscripts,  as  for  example 


3    8"^w  .  -  -  ...  .,.,  ^. 


3x3  y 


The  volume  element  for  physical  space  and  parameter  space  will  be  respectively  dX 
and  dU.  We  define  a  set  of  rules  not  unlike  summation  notation  to  shorten 
expressions.  First,  stand  alone  variables  or  powers  of  stand  alone  variables  are 
summed  over  all  free  indices.   As  an  example 

2    3   3      2 

u^      =     I        I         u.^         .  (III. 6) 

i=l  j=l 

For  products  of  two  or  more  variables  summation  is  done  over  common  Indices.  For 
example 

3   3 

i=l  j=l 

Parentheses  are  often  used  to  indicate  that  summation  is  done  first  before  any 
other  operation.   For  example 

(uj^Ujb)2  =  (  ^  u/ujb)2   .  (III.,:) 

J 

With  the  notation  defined,  the  variational  principle  for  the  mesh  generator 
can  be  defined.  The  smoothing  functional  extends  the  two  dimensional  case  by 
adding  the  magnitude  of  the  gradient  squared  of  the  third  parameter  coordinate. 


^s  ~  /  "j^   dX   .  (III. 9) 


The  weighting  functional  has  exactly  the  same  form  as  In  two  dimensions 
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I„  =  /  wJ  dX 


(III. 10) 


where  J  Is  the  Jacoblan  In  three  dimensions 


J  =  det 


^u   ^   ^. 


Y    Y    Y 
u    V    w 


Zu    ^    ^w 


(III. 11) 


The  Jacoblan  can  be  expressed  In  terms  of  the  Indexed  notation 


■^   ~   ^Ijk'^i  '^j  '^k 


(III. 12) 


where  e  ^ -ji^^  is  the  permutation  symbol.  The  orthogonality  functional  is  more 
complicated  In  three  dimensions  because  it  must  measure  the  orthogonality  of 
three  sets  of  grid  lines  Instead  of  one 


1„  2x2  ^  ,.,  1„  3s2  ^  ,.,   2.,  3x2,  ,3 


I^  =  /  [(uj'uj^)^  +  (Uj^uj^)^  +  (Uj^u.-^)^] 


J-'dx 


(III. 13) 


Next,  as  before,   the  variational  principle  can  be  written  in  terms  of  the 
parameter  coordinate  u,v,w 


Ig  =  /  u^      Jdu 


(III. 14) 


1  ,  =  /  wJ^  du 
w   ^ 


(III. 15) 


1..  2x2  J.  /..  1.  3x2  ,  /..  2..  3x2,  tA 


lo  =  /  [("j  "j^)  +  ^"J  "j^  ■"  ^"j  "j^  5 


J^  du   . 


(III. 16) 
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Determlnatlon  of  the  Partial  Differential  Equations 

With  the  Integrands  written  with  the  paramater  variables  as  the  independent 
variables  the  three  Euler  equations  can  be  written  down.  Using  the  summation 
notation  the  Euler  operator  is  (as  before,  summation  is  performed  on  the  matched 
Indices  k  and  1  multiplies  9 /3 x^  to  prevent  summation  on  £ ) 


Lf  =  1  •  -L-  -  -i ^f  =  0   il  =  1,  2,  3  (III. 17) 

a/    3u^  3x^^ 


where  f  is  the  weighted  sum  of  the  integrands 

f  =  Xg(u/  )J  +  X^wj2  +  X^[(uj^uj2)2  +  (^,_^1^,^3)2  +  (^^2^^3)2j  jA   ^  (m.ig) 

To  facilitate  calculation  of  the  equation  several  Lemmas  will  be  introduced. 
Lemma  1 


Ju-   =  


8X3 


Lemma  2 


T   i    1-9-^ 
J^j^ik^  =  -17 


Lemma  3 


3u^  8x^J 
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Lerama  4 


3    3 
J  =  — —  3 


3u^  3x^j     3x^J  3u^ 


Lemma  1  follows  from  the  chain  rule,  Cramer's  rule  and  the  rule  for 
differentiating  a  determinant.  It  gives  the  explicit  representation  of 
u^   in  terms  of  functions  with  dependent  variables  u  . 

Lemma  2  follows  from  the  chain  rule  and  the  rule  for  differentiation  a 
determinant. 

Lemma  3  follows  from  direct  calculation  using  the  summation  shorthand  for  J. 

Lemma  4  as  does  Lemma  3  follows  from  direct  calculation. 


Derivation  of  Smoothing  Functional 

To  start,  notice  that  u^  J  does  not  explicity  contain  x*^  so  that 


a      2 
-^  (u.^  J)  =  0   .  (III. 19) 

3/ 


Hence  the  Euler  equation  becomes 


^    ^   (u/  J)  =  0  (III. 20) 


0  =  - 


9  fo..  1  ^""Z  T  ^ ...  i^  3 J 


(2u/  —2^  J  +  (uji  )  ^)   .  (III. 21) 


3u^    ^  3x1^^       ^    3X1^^ 
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(using  Lemma  1) 

,2 


(again  using  Lemma  1  and  simplifying) 


3u^      ^    9x1,^     ^  3x1,^  Bx^J 


(differentiating,  using  Lemma  3,  and  collecting  terms) 


^  J  ^     9X1^^  3xiJ    3u^ 


By  using  Lemmas  1  and  2 


9"j^       1„   n„  m  ,  1   9    9 ^   .  (TII.2S) 


3 „  i^  %  ™  + 


Combining  this  with  the  previous  expression 


„.-.(„/KVv."-^-/^."''n"^^ 


.   I,  8   3J    2  3    3J    3    3J 

+  2u.^uj,  "^  — ; 7  -  ^  — 7-   7  r  7' 

^  3u^  3xiJ   ^  ^u^   3xiJ  3x/  Sx^J 


(III. 26) 
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(uslng  Lemma  2  and  combining  terms) 


8x1^  3x^J 


-U,  Uo    r   +  -   r  —0 r  J 


'j  "ii 


ax/3x,J   Jax„"8x,J8x^^3x,J 


Derivation  of  Weighting  Functional 

As  before,  the  Euler  equation  is  written, 


3/    3u^  3xj^^ 


Carrying  out  the  differentiation  and  using  the  chain  rule  yields 


n    t2  3w    k   OT   9    9  J    o   3J   3J      ..  3J   3w 
0  =  J  —   U£   -  2Jw  — ; _  -  2w  ; —   J  -  2J 


k  2..k  a„^£       9^^£  g^k        g^^£  g^^k 


(III. 27) 


0  =  [  1  •  J_  -  J —   )  wj2   £  =  1,  2,  3  (III. 28) 


Finally,  using  Lemmas  3,  2,  1,  and  simplifying  gives 


0=  -2\k^A  V  -  A^  ^  •  ("1-29) 

3u'^ 
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Derlvatlon  of  Orthogonality  Functional 

It  suffices  for  our  purposes  to  calculate  the  Euler  Equations  for  the 
functional 

/  (Uj3u>)2j^  du  .  (Til. 30) 


Again  since  the  integrand  does  not  depend  explicitly  on  x  ,  the  Euler  equation 


can  be  written, 


^    ^   (u,^u  b)2jA  (III. 31) 


au'^axk^  '    ' 


Using  Lemma  1  yields 


2 

0  =  -  ^_i_(    Ji__li_)  (III. 32) 

3u"  ax/        Sxgj  3xijJ 


Differentiating   gives 


0   =   -   2^[         '^'        J^l^^        ^'J        IL.1^1^].  (III. 33) 

3uk     sx^'^xj^^  ax^J  ax^J     3xt,'=ax/  3xJ  ax^J  ax^™ 


Then,    by  Lemma   3, 


0  =  -   9^J   a  r     aj  3j  3J  >,  _   3^j   a  ,     aj  aj  aj  ^      ^^ 
ax^'^ax,^  au'^  3x  ™  ax^j  ax,  J    ax.'^x  ^  au'"  ax  J  ax,  J  ax  " 

aK  Dab        DU  aba 


By  using  Lemmas  4,  2,  and  then  1,  we  find, 
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3        3J             9        8J            9              i      1,               i      9         9J  ,^^^    ,,. 

— i = 57  °  '^ik  "i^  "  ^ik r     •  (III. 35) 


Finally  by  using  the   two  end  expressions   In  the  previous  equation 


29  2j        rf9J       SJ^  n        3^J        U.3J       3J  n        3^J 


0  =   --^^^[(    -^^)    X,," 


Xj 


'-aW    ^^    '-J'-^^'      '''     9-b"^-l"       S'^b'^^V      ''^     S-a^^V 


Sxb'^^^a^  3xbJ3xi" 


(III. 36) 


U      7    tJ       Xj^"    _ -    + 


Ik    : — - — r  +  : — r  : — i  Hk 


Sxb'^x^^  ax^JSx^J  axa'^xi"        3x3">3x3J  9x^,^9x1" 


3x3-3x,J      ''     Sx^JSx," 


(To    be    consistent   with    the   previous    notation  m   will    be    exhanged    with    1    and    then 
Lemma   1    will   be  used   to   simplify) 
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^       b     b        3^-J         .        b     a        8^J 

am  Dm 


(III. 37) 


,         SJ         fab         S-'  ,aa         9J 


3x,i3x/        ^      ^     3x,i3x^"  '     3x,J3x/ 


By  setting  (a,b)  =  (1,2),  (1,3),  (2,3)  and  summing  the  expressions  the 
coefficients  of  the  orthogonality  contribution  to  the  PDE  system  can  he 
evaluated. 

The   form  of    the   coefficients    is   composed   of   only   two   building   blocks 


i        ^         9    -J 


'  3.i«x/ 


Uj  ^  is  a  set  of  nine  expressions  while  3  J/3x.J3x,  is  a  set  of  81  expressions. 
However,  3  J/3x.-J9x,  has  only  18  unique  nonzero-terms  or  36  nonunique  nonzero 
terras. 

By  summing  on  i  and  j  the  coefficients  of  the  second  order  differential 
operator  terms  can  be  calculated.  However,  it  is  not  necessary  even  to  evaluate 
all  the  coefficients  since  another  symmetry  exists.  By  exchanging  £  with  n  and  k 
with  m  it  is  seen  that  the  system  of  equations  is  symmetric  and  so  only  the  lower 
or  upper  triangle  of  coefficients  of  the  second  derivatives  need  by  evaluated. 
That  is  the  coefficients  can  either  be  written 
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^11  ^^12    •••   ^^32  ^33  yil   yi2    •••   y32  ^33  ^11   ^12    ••'    ^32   ^33 

yii  yi2  •••  y32  ^33  ^ii  ^12  •••  ^32  ^33 

^11  ^12  •••  ^32  ^33 
or 

'^ll  ^^12  •  ••  ^32  ^^33 

^11  ^12  •••  ^32  ^^33       yil  yl2  '••  ^32  ^33 

x^l  xi2  •••  ^^32  ^^33       yn  yl2  •••  ^32  y33        ^11  ^12  '••  ^32  ^33 

There  is  still  the  question  of  normalization.   It  one  were  to  choose  X  ,  X  ,  and 

X^  all  of  the  order  of  one  the  smoothing  functional  would  dominate  the  solution 

algorithm.   That  is  to  say  changes  in  the  mesh  spacing  fron  the  orthogonality 

coefficients  or  weighting  coefficients  would  be  small  compared  to  those  of  the 

smoothing   coefficients.    This   was   observed   in   performing   the   numerical 

calculations . 

It  is  desireable  to  have  a  normalization  (a  constant  multiplier  for  each  X) 

so  that  noticeable  changes  in  the  mesh  spacing  are  caused  by  the  coefficients  of 

the  corresponding  functional  when  the  associated  X  is  of  order  one.   As  before, 

in  two  dimensions,  using  simple  dimensional  analysis,  constant  multipliers  X  ', 

X,,',  and  X^'  can  be  found  that  seem  to  give  the  desired  normalization.   This 
wo 

gives  us 


^.  .    .2      X 

I  =   ^ 


^s' 


/  u^J  dx  +  ^  /  [(Ujlu.2)2  +  (Ujlu.3)2  +  (..2^.3)2,  j3  ,^ 
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(III. 38) 


+  — ^  /  wJ  dX   . 


^o 


X  '  is  set  to  one  and  X^'  and  X^'  are  determined  so  that  all  three  functionals 
have  the  same  dimensions.  The  integrand  of  the  smoothing  functional  has  units 
parameter  space  length  I    over  physical  space  length  £  squared 


I^~  (  ^)  (III. 39) 


while 


^.~r,f 


I„~  w[  i)^  (III.AO) 


where  w  are  units  of  the  weight  function  w.  This  implies 


X^'~  [  i)^  (III. 41) 


^v'~  [^  t  -)  ]  (III. 42) 
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so  that  all  the  units  will  agree. 

The  scale  length  £  is  chosen  to  be  the  average  length  scale  of  the  physical 
space  while  £  is  the  average  of  the  number  of  mesh  points  in  each  of  the 
parameters  space  coordinates,  w  is  chosen  to  be  some  approximation  to  the  volume 
average  of  the  weight  function  over  the  physical  domain.   That  is 


w  =  i  /  wdX  (III. 43) 

V 


such  that 

V  =  /  dX  .,,.,,    ,  .  (III. 44) 

where  ^    is  the  physical  domain. 

With  the  system  of  equations  0=1    a.  .  x.  .   +  i,    written  down,  it  is  now 

i.j.k   iJ   iJ     b 

possible  to  discuss  the  iterative  numerical  solution  of  these  equations.  The 
first  step  in  the  discussion  is  to  define  the  physical  and  parameter  domains. 
For  now  the  physical  space  will  be  a  simply  connected  region.  The  parameter 
space  will  be  a  cube  in  Eucidean  3-space  with  sides  of  lengths  IMAX-1 ,  JMAX-1 , 
KMAX-1.  IMAX,  JMAX,  and  KMAX  are  respectively  the  number  of  points  in  the  u,  v, 
and  w  directions  used  in  solving  the  finite  difference  equations. 

With  the  parameter  space  being  a  cube,  all  differencing  of  the  first  and 
second  derivatives  can  be  done  using  centered  differences.  For  example,  at  the 
point  (u,v,w)  =  (i,j,k)  and  iteration  n 


c(i  +  l,j,l 


3x  _  x(i  +  l,j,k)  -  x(i  -  l,j,k) 

5—-  5 (III. 45) 

9  u  2 
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s2„     (n)  (n)        (n) 

L^  =  y(i.j.k+l)  -  2y(i.j.k)  +  y(l.j,k-l)  (III. 46) 


3w2 


,      (n)  (n)  Cn)  fn) 

l\_  ^   z(l  +  l.j+l,k)  +  z(i-l.j-l.k)  -  z(l  +  l.j-l.k)  -  z(l-l,j+l.k)   ^  (III. 47) 


8  u3  V 


As  before,  the  residuals  of  the  three  equations  are  written  down 

^rCti)  =   \  ^^ii'^'^ij^  "^  ^^   £  =  1,  2,  3  (III. 48) 

i,j,k 

where   a^^^^  or  the  coefficients  of  the  2nd  derivatives  calculated  proviously  with 

n 

the  derivatives  replaced  by  the  differences.   The   b  are  terms  without  2nd 

derivatives  and  are  calculated  using  differences  in  place  of  derivatives.   Thus, 

at  every  interior  point  of  the  cube  the  following  system  is  solved  for  the 
(n+l)th  iterate. 

o  /  X    r.         J    J  (n+1)  J  (n) 

0  =  ^R(n)  +  I      -  2ajjMxi     (i,j,k)  -  xi    (i,j,k)]   .         (III. 49) 

i.j 

If  x"(i,j,k)  is  updated  iminediately  so  that  the  new  value  is  used  in  the 
iteration,  then  the  method  of  iteration  will  be  called  Gauss-Seidel.  If 
x"(i,j,k)  is  not  updated  the  iteration  is  called  a  Jacobi  iteration. 


Computational  Examples 

In  the  following  sections  sample  computations  will  be  given  demonstrating 
various  characteristics  of  the  variational  principle  and  its  solution  algorithm. 
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Welghtlng  Functional  Test 

As  a  first  example,  the  parameter  space  cube  will  be  mapped  one  to  one  onto 
a  unit  cube  at  the  origin  In  physical  space  (Fig.  25).  The  initial  guess  or  the 
numerical  iteration  will  be  a  lattice  of  points  uniformly  spaced  in  all  three 
directions.  The  computation  uses  20  mesh  points  in  each  parameter  space 
direction  (IMAX  =  JMAX  =  KMAX  =  20)  for  a  total  of  8000  points.  However  only 
(18)  =  5832  points  are  iterated  upon  since  the  other  points  are  fixed  at  the 
boundary. 

If  X  =  1  and  X^  =  Xq  =  0  it  is  seen  either  from  direct  computation  or  from 
the  continuous  equations  that  the  initial  guess  is  a  solution  of  both  the 
continuous  and  discrete  case.  However,  if  one  specifies  a  weight  function  and 
Increases  X  away  from  zero  the  initial  guess  is  no  longer  a  solution.  As  an 
example,  the  effect  of  the  following  weight  function  will  be  examined, 

w(x,y,z)  =  1000  exp{(0.25  -  (x^  +  y^  +  2^)1/2)2/0. 05]   .  (III. 50) 

Clearly  the  effect  of  the  weight  function  should  be  to  bunch  the  mesh  around  a 
spherical  shell  centered  at  the  point  (x,y,z)  =  (0.5,0.5,0.5).  Table  5  compares 
the  numerical  value  of  the  smoothness  functional  and  weighting  functional  as  X^ 
varies  between  0  and  1  for  a  fixed  number  of  iterations.  This  shows  that  the 
weighting  functional  is  indeed  minimized  as  X  is  increased.  The  graph  In 
Fig.  26  summarizes  the  results  of  the  table.  To  illustrate  how  the  mesh  has  been 
deformed  various  planar  cuts  have  been  made  in  the  physical  cube  (Figs.  27a, b). 
All  the  cuts  hold  one  parameter  variable  fixed  and  vary  the  other  two.  On  each 
cut  one  can  see  how  the  lines  are  deformed  on  the  plane  In  a  circular  manner  and 
how  the  plane  itself  is  deformed  in  the  way  one  would  anticipate. 


-66- 

( TABLE  5) 

Number  of  iterations:  50 

X  Weighting  Functional      Smoothness  Functional 

(Normalized) 

0.0001           3.045  X  10^  9.208  x  10^ 

0.000316          2.979  x  10^  9.207  x  10^ 

0.001             2.820  X  10^  9.208  x  10- 

0.00316          2.539  x  10^  9.227  x  10^ 

0.01             2.157  X  10'^  9.301  x  10^ 


0.0316  1.815  X  10^  9. 442  x  10^ 

1.585  X  10^  9.642  x  10^ 

1.473  X  10^  9.880  x  10^ 

1.439  X  10^  1.011  X  10^ 


1.0 

Orthogonality  Test  of  Functional 

This  test  problem  will  map  one  to  one  the  parameter  cube  onto  a  frustum  as 
Illustrated  in  Fig.  28.  Apaln,  twenty  mesh  points  will  be  used  in  each  direction 
in  parameter  space. 

The  coordinates  of  the  bottom  plane  of  the  frustrum  are 
(x,y,2)  =  (0.25,0.25,0.5),  (0.25,0.75,0.5),  (0.75,0.75,0.5),  (0.75,0.25,0.5) 
while  the  coordinates  of  the  top  plane  are  (x,y,z)  =  (0,0,1),  (1,0,1),  (1,1,1), 
(0,1,1). 

To  generate  an  initial  mesh,  the  following  procedure  is  used 

a)  Twenty  points  are  uniformly  distributed  on  each  edge  of  the  frustrum  as  in 
Fig.  29  (each  corner  of  the  figure  should  contain  a  mesh  point). 

b)  On  each  face,  the  points  on  opposite  edges  are  connected  to  give  a  grid  on 
each  face  (Fig.  30).  (It  is  understood  that  parallel  lines  in  parameter 
space  do  not  intersect  in  physical  space.   edges  should  not  intersect.) 

c)  The  points  of  intersection  on  opposite  faces  are  connected  to  give  a  3-D 
lattice.  (It  is  understood  that  lines  emanating  and  landing  from  and  onto 
the  same  face  should  only  form  rectangular  columns.) 

The  last  step  is  hard  to  illustrate  but  hopefully  easily  understood. 
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As  with  the  weight  functional  example,  if  ^q  =  ^u  =  0  a"<^  ^s  ~  ^  both  the 
discrete  and  continuous  sets  of  equations  are  solved  with  the  Initial  mesh.  By 
increasing  X^  however,  the  solution  of  the  finite  difference  scheme  starts  to 
change.  Table  6  compares  the  numerical  values  of  the  smoothness  functional  and 
the  weighting  functional  as  X^  varies  between  0.1  and  10  for  a  fixed  number  of 
iterations.  Figure  31  summarizes  the  results  of  table  6.  Once  again  the 
orthogonality  functional  is  minimized  at  the  cost  of  some  decrease  in  smoothness. 
Two  planar  cuts  are  also  (Fig.  32)  given  below  to  demonstrate  the  nature  of  the 
solution  as  orthogonality  weighting  is  increased.  Each  horizonal  cut 
(w  =  constant)  shows  a  bending  of  the  whole  plane  to  minimize  the  orthogonality 
integral. 

A  nontrivial  example:   Zoning  a  section  of  an  aircraft  fuselage  and  wing 

One  area  in  which  3-D  meshes  are  used  is  in  mathematical  modeling  of  air 
flow  around  an  aircraft  fuselage.  A  fuselage-wing  section  is  illustrated  in 
Fig.  33.         . 


(TABLE  6) 

Number  of  iterations; 

Orthogonality  Functional 
(Normalized) 


50 


Smoothness  Functional 


0.1 

0.3363  X 

10^ 

0.6844  X  10^ 

0.316 

0.3319  X 

10^ 

0.6848  X  10^ 

1.0 

0.3197  X 

10^ 

0.6858  X  10^ 

3.16 

0.2924  X 

10^ 

0.6881  X  10^ 

10.0 

0.2A60  X 

10^ 

0.6925  X  10^ 

31.6 

0.1902  X 

10^ 

0.6992  X  10^ 

100.0 

0.1456  X 

10^ 

0.7074  X  10^ 

316.0 

0.1212  X 

10^ 

0.7158  X  10^ 

1000.0 

0.1110  X 

10^ 

0.7232  X  10^ 
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As  can  be  seen  from  the  complexity  of  the  geometry,  a  detailed  description 
of  the  location  of  surfaces  is  not  helpful.  However  a  discussion  of  the  general 
properties  of  the  mesh  would  be  helpful.  To  understand  the  logical  arrangeraent 
of  the  parts  of  the  mesh  it  is  best  to  start  with  the  unit  cube  in  parameter 
space  mapped  onto  a  parallelepiped  and  then  deform  the  parallelepiped  in  steps 
(Fig.  34).  For  clarity,  hidden  sides  are  shown  with  dotted  lines  and  each  corner 
of  the  figures  is  numbered  in  an  appropriate  fashion.  The  first  step  is  to 
deforn  the  parallelepiped  into  a  "c"  as  shown  in  Fig.  35.  Next,  the  surfaces 
with  corner  at  1,A,8,5  and  2,3,6,7  are  bent  to  conform  to  the  fuselage  and 
"glued"  to  the  fuselage  (Fig.  36).  Finally  the  c  is  pushed  together  around  the 
wing  as  shown  in  Fig.  37.  The  number  of  mesh  points  in  the  u  direction  is  chosen 
to  be  twenty-three  (imax  =  23).  The  number  of  mesh  points  in  the  v  and  w 
directions  is  chosen  to  be  ten  (jmax  =  kmax  =  10). 

To  generate  an  interior  lattice  a  procedure  similar  to  the  one  used  by  the 
orthogonality  test  problem  is  used.  That  is,  boundary  lines  are  subdivided,  then 
boundary  surfaces  are  generated,  and  finally  a  lattice  is  constructed  from  the 
grid  surfaces.  Several  cuts  of  the  lattice  are  displayed  in  Fig.  38  to 
illustrate  the  initial  setup. 

With  an  initial  guess  generated  the  smoothing  parameter  X^  is  set  to  one  and 
X  =  X  =  0.  Figure  39  shows  the  solution  after  50  iterations.  The  most 
Interesting  feature  of  the  solution  is  the  way  the  mesh  is  concentrated  around 
the  wingtip.  This  can  be  explained  by  fact  that  three  potential  equations 
(Laplace)  are  being  solved.  By  an  analogy  with  electrostatics,  one  can  see  that 
the  wing  is  acting  as  a  lightning  rod  in  the  x(u)  direction  with  a  corresponding 
concentration  of  level  curves.  That  is,  equipotentials  become  close  together 
near  the  sharp  point. 
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For  many  problems,  the  grid  as  it  Is,  is  satisfactory.  However,  there  might 
be  other  features  of  the  flow  that  need  to  be  resolved  as  well.  To  demonstrate 
the  capabilities  of  the  weight  functional,  to  alter  the  grid  as  needed,  a 
weighting  which  is  large  in  a  region  away  from  the  wing  tip  will  be  used.  To  be 
more  precise,  the  weight  function  will  be  large  in  a  an  elliptic  region  far  from 
the  tip  of  the  wing  as  shown  in  Fig.  AO. 

The  region  is  removed  far  enough  from  the  wing  tip  to  pull  the  mesh  away 
from  the  point.  The  set  of  cross  sections  in  Fig.  Al  illustrate  the  effect  of 
the  weight  function. 
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IV.  Applications  in  Two  Dimensions 

In  this  chapter  two  problems  are  solved  numerically.  They  are  the  evolution 
of  the  Raylelgh-Taylor  instability  and  the  step  in  a  Mach  3.0  wind  tunnel 
problem.  Specific  regimes  of  these  problem  will  discussed  later.  Immediately 
below  the  equations  governing  these  problems  and  approximation  of  these  equations 
are  discussed. 

The  equations  common  to  the  two  problems  are  the  inviscid  compressible  flow 
equations.   They  are,  using  a  Lagrangean  formulation, 


^  +  pV  .  u  =  0  (IV.  1) 

dt 


p  ^  +  Vp  =  0  (r'.2) 

d  t 


p  ^  +  pV  .  u  =  0  (IV .  3 ) 

dt    '^ 


p  =  (Y  -  Dpi   .  (IV. A) 

The  notation  is  standard,  but  for  completeness  the  definition  of  each  s>Tnhol  is 
given  below.   For  more  information  about  these  equations  see  [5B,  7B,  and  8B] . 

symbol                 definition  units 

p           density  (m/£  ) 

t           t Ime  ( t ) 

u           fluid  velocity  vector  (Jl/t) 

p          pressure  (m/£t  ) 
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I  internal  energy  per  unit  mass         (&  /t  ) 

Y  ratio  of  specific  heats  (unltless) 

(m  =  mass,  £  =  length,  t  =  time) 

To  numerically  represent  and  solve  these  equations  a  von  Neumann-Rich tmyer 
type  formulation  is  used  [IB].  That  is,  the  variables  are  staggered  on  the 
computation  mesh.  The  velocities  are  calculated  at  vertices  formed  by 
intersection  of  mesh  lines  and  are  called  vertex  centered  quantities.  The 
internal  energy  and  density  are  calculated  at  the  center  of  a  four  sided  polygon 
formed  by  four  mesh  lines.  This  polygon  will  be  called  a  cell  and  the  internal 
energy  and  density  are  called  cell  centered  quantities. 

Representation  of  cell  and  vertex  centered  quantities  are  needed  everywhere. 
To  do  so  cell  centered  quanties  will  be  assumed  to  be  piecewise  constant  In  each 
cell  and  vertex  centered  quantities  are  represented  using  a  piecewise  bilinear 
approximation.  For  each  cell  there  are  four  vertices  as  in  Fig.  A2.  The  cell 
can  be  mapped  to  a  unit  square  as  shown  in  the  same  figure  using  the  mapping 

X  =  C(l  -  n)x^  +  Cnx2  +  (1  -  nnx3  +  (1  -  C)(l-n)x4  (TV. 5) 

y  =  C(l  -  n)yi  +  Cny2  +  (1  -  C)ny3  +  (1  -  ?)(1  -  n)y4  (IV. 6) 

4.  =  5(1  -  n)(j>i  +  5n<t>2  +  (1  -  C)n(t>3  +  (1  -  C)(l  -  n)4'4  .         (iv.7) 

This  mapping  is  convenient  in  calculating  many  things  [ 14A] .   To  calculate  a  cell 
area  the  integration  has  been  performed  on  the  unit  square. 
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V,  "jj    dxdy  =/   /  l^i^dl,dy      .  (IV. 8) 

The  cell  area  is  found  to  be  the  sum  of  four  subareas  associated  with  each  vertex 
of  the  cell 

-c   =  c^l   ^  ch  -^  ch  ^  c^  ^^^'-^^ 


where 


^Ai    =  -i   [(X2   -   xi)(y^   -  Yi)    -    (72   "  YiXx^   -  x^)]  (IV. 10) 


c^2   =  -J-   [^^^3   "  '^2)(yi   ■  ^2^   ~   (^3   "  y2>('^l   "  ""2^^  (IV. 11) 


ch   =  i    ^('^A    -   ^3^^y2    -   ^3^    -    ^^4    -   y3)<^2   "   ^3^1  (IV. 12) 


c^A   =  7^   t(xi   -  X4)(y3   -  y^)   -   (y^   -  y4)(x3   -  x^)]  (IV. 13) 


Discretlzing  the  Equations  of  Motion 

In  discretlzing  (IV. 1)  to  (IV. A)  derivatives  of  vertex  centered  quantities 
at  cell  centers  and  derivatives  of  cell  centered  quantities  at  vertices  are 
needed.  Cell  centered  derivatives  at  vertex  quantities  are  calculated  using  the 
average  value  over  the  entire  cell 
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<  I7  >  ■  ^ // 17  "-"^  •  <"-^^> 


As  with  the  subareas  the  approximate  derivatives  can  be  represented  as  a  sun  of 
terms  associated  with  the  vertices  of  the  cell.  These  coefficients  will  be 
called  geometric  coefficients. 

^  9^  ^  "  \r  U'^'^l*!  +  cC'^2*2  +  c'^^S^B  +  cC'^A^J  (IV. 16) 

c 

^  T^  ^  '  v~  ^c^yi*i  "^  ccy2*2  +  c<^y3*3  +  c^y^^J  (iv.i?) 

"      c 


ccxi  =  _  (72  -  74)  (IV. 18) 


cCX2  =  2  <^y3  "  yi)  (IV. 19) 


1  / 

cCX3  =  -  (y^  -  72)  (IV. 20) 


c^H   =  7  (yi  "  ^3^  (IV. 21) 


and 


-74- 


ccyi  =  i  (X4  -  X2)  (IV. 22) 


c 


cy2=i(xi-X3)  (IV. 23) 


,cy3  =i  (X2  -  x^)  (IV. 24) 


ecy4  =^  (X3  -  xi)  (IV. 25) 

Again  notice  the  cyclic  pennutatlon  of  the  indices.  As  it  turns  out  these 
coefficients  can  be  used  to  calculate  the  derivatives  of  cell  centered  quantities 
at  vertices.  To  generate  derivatives  of  cell  centered  quantities  at  vertex 
centers  we  devise  a  quadrature  formula  using  the  subareas  ^A^ ,  ^^2^  c^3 '  ^"^  c^4 
[7A].   To  integrate  a  vertex  centered  quantity  over  a  cell  the  approximation 

/   (t)  dx  =  ^A^^^   +  ^A24)2  +  0^3*3  +  cMa        ''  (IV.  26) 

cell 

will  be  used.  With  (IV. 16),  (IV. 17),  and  (IV. 26)  geometric  coefficients  for 
vertex  centered  can  be  found.  Consider  the  following  set  of  four  mesh  cells  in 
Fig.  43.  Let  the  p^  be  cell  centered  quantities  and  let  u  =  (u,v)  be  a  vertex 
centered  vector  quantity.  Also  assume  at  the  exterior  vertices  the  vertex 
centered  vector  quantities  are  zero.  Then  since  u  is  zero  identicallv  on  the 
boundary  of  the  four  cells,  it  follows  that 

/  7  •  (pti)  dx  =  ^  pti  •  dn  =  0   .  (IV. 27) 
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Expanding  (IV. 27),  we  write 

/  (pV  •  u  +  u  •  7p)  dx  =  0   .  (IV. 28) 

Using  (IV. 26),  (IV. 16),  (IV. 17)  and  assuming  Vp  is  a  vertex  centered  quantity  in 
IV. 28,  we  obtain 

0  =  ^IPI  [  T5-  (1CX3U  +  icy3v)]  +  V2P2  [  y-  (2cx^u  +  2cy4v)] 
^1  ^2     ^ 

■*"  ^'3P3  [  y-  (3'=^1^  "^  3<=yi^)]  "^  ^'4P4  [  y-  (^^^2"  ^   4^y2^)] 


4 


+  "  1^  (1A3  +  1A4  +  lAi  +  1A2) 


9x 


dp 


(IV. 29) 


+  ^  9^  (1A3  +  iAa  +  lAl  +  1A2) 


If  u  =  1  and  V  =  0  then  3p/3x  is  isolated  as 


3p  _  -Pi  1^^3  "  P2  2^x4  -  P3  3^^1  -  P4  4^^2  (IV  30) 

3x  jAj^  +  ]^A2  +  1A3  +  ^A^ 


Similarly  if  u  =  0  and  v  =  1  then  we  can  write 


3p  _  -Pi  l^y3  -  P2  2'^y4  -  P3  3^yi  -  P4  4^y2  (IV  31) 


3y  iM  "^  1^2  "^  1^3  "•"  1^4 
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By  letting 

^v  "   1^1  +  1^2  +  1^3  +  1A4   .  (IV. 32) 

We  find  that  differentiation  formulas  have  the  same  form  as  (IV. 16)  and  (IV. 17). 

With  the  above  apparatus  in  place  it  is  now  possible  to  discretize  (IV. 1)  to 
(IV. 4).  Vertex  quantities  are  integer  indexed.  For  example  Uj ^  is  the  x 
component  of  the  velocity  field.  Cell  centered  quantities  will  be  half  integer 
indexed.  For  example,  V±+if2  1+1/2  ^^  ^^^  pressure  at  the  cell  center  (i  +  1/2, 
j  +  1/2).  Left  superscripts  indicate  the  time  level.  "'^"i-i"  Is  the  x  component 
of  the  velocity  field  at  the  n   time  level.  At  is  the  time  increment. 

At  time  level  n  the  following  quantities  "1^+1/2  i+l/2»  ^^i+1/2  i+1/2'  '^'^ii' 
"v..,  and  '^m.  .  are  stored  where  ^^i    is  a  vertex  mass.   The  first  step  in  the 


calculation  is  to  find  the  pressures  using  (IV. A) 

"Pi+l/2,j+l/2  =  (Y  -  1)  Pi+l/2,j+l/2  Ii+l/2,j+l/2  +  Pi+l/2,j+l/2   '  (1^.33) 

Pi+1/2  i+1/2  ^^    ^    viscous  pressure  to  be  defined  later.   Next,  the  fluid  is 
accelerated  using  the  pressures 


"m.  .  (""^^u.  .  -  "u.  .)     .n 

_^^^ ^i ii:  +  fL_P)    =0  (IV. 34) 

%,        At  ax  \. 


where  "v. .  is  a  vertex  volume  as  defined  in  (IV. 32)  and  the  derivative  of  p  is 
found  using  (IV. 30).   Similarly 
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%j  ("^Sj  -  Nj)  ^ .  3- 


"V,  ^  At  ^^       ij 


[    t^]        -  0      .  (IV. 35) 


ij 


Next,    the   internal  energy  is  updated  using  the   following  discretization. 


/n+lf  _  Ht  \ 

n+1  ^        H-H/2.J+1/2  ^i-H/2.j-H/2^ 

Pi+l/2,j+l/2  n (I^.36) 


Pi+l/2,j+l/2  ^        "i+l/2,j+l/2)  -  ° 


where 


ij   2  ^   iJ       iJ 


The  fluid  mass  is  a  Lagrangean  invariant  and  so 


n+^m. j  =  "m. j   .  (IV. 37) 


Next,   the   georaetry   is   updated,   for   a   Lagrangean   calculation,   using   the 
approximation 


At  ij 


With  (IV. 39)  the  geometric  coefficients  and  volumes  may  be  found. 


Finally  "  ^1+1/2  i+1/2  ^^    found  using 
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n+1  f       Ij  A     J-     i+lJ 

Pi+l/2,j+l/2   =  I    vT:  i+l/2,j+l/2A4  +  v..,     .   i+1/2 , J+I/2A1 

1 J  1+1 , J 

(IV. 39) 
•^  V. ,T     ..,    1+1/2, J+I/2A2  +  V.     ., 1    1+1/2,^+1/2^3   J/^i+l/2,j+l/2 

l"ri  J  J"r  i  1  ,  J  ""i- 

At  this  point  several  comments  are  appropriate.  First,  boundary  conditions  are 
easily  implemented.  Free  slip  boundaries  only  require  that  the  acceleration  of 
the  vertices  be  tangent  to  a  boundary.  Secondly  mass,  momentum,  and  energy  are 
all  conserved.  Mass  is  trivially  conserved  because  of  (IV. 37).  Momentum  is 
conserved  because  of  the  way  the  geometric  coefficients  were  defined  in  (IV. 30) 
and  (IV. 31)  via  (IV. 28)  [lAA].  Energy  conservation  is  proven  as  in  the 
continuous  case  by  dotting  "  '  u .  .  into  (IV. 34)  and  (IV. 35),  then  adding  the 
result  to  (IV. 36). 

The  viscous  pressure  term  in  (IV. 34)  is  of  the  form 


*                 1  ^i+l/2,j+l/2   .   ,n  ^    -►^  ^  /^  ^2  frv   /n^ 

Pi+l/2,j+l/2  =  "  ^  A  -KT """^  ^  '^  *  "^i+l/2,j+l/2)  (^'^^   (IV. 40) 


2 
where  w  is  a  positive  number  between  zero  and  one  and  (Ax)   has  units  of  length 

2 
squared.   (Ax)   is  taken  to  be 


(Ax)2  =  [  [   _^J^ '^      ^        ]  .  (IV. 41) 

1=1        (V^2) 
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A  subsequent  remap  phase  can  take  the  data  at  time  level  n+1  and  map  it  onto 
a  different  grid  in  which  case  (IV. 38)  will  be  replaced  by  another  rule.  Because 
of  the  relative  motion  of  the  mesh  and  fluid,  the  fluid  streams  through  the 
computation  mesh.  The  governing  equation  for  this  streaming  is  the  following 
integral  equation, 


__//(})  dxdy  +  ^  (J>u  •  dn  =  0   .  (IV. 42) 


The  equation  states  that  a  conserved  quantity  (})  in  a  region  n  changes  only  from 
contributions  flowing  in  from  the  boundary  3f2 .  The  inflow  is  controlled  by  the 
vector  field  u.   If  $>  is  dif  f  erentiable ,  the  divergence  theorem  yields 


^  +  V  •  (4.U)  =  0   .  (IV. A4) 

o  t 


To  approximate  the  flux  term  of  (IV. A3)  flux  velocities  are  calculated  at 
all  the  cell  vertices. 


."n  =  "^^""ij  -  ^^^^T^  -  ""'^.  -  ."ij   • 


Notice  that  if  the  mesh  moves  with  the  lagrangean  velocity  x.  •  =  "''"■'^x.  .,  then 
there  is  no  flux.   This  means  no  remap  is  necessary. 

To  determine  the  new  mesh  the  grid  generator  described  in  chapter  II  is  used 

in  the  same  manner  as  was  the  one-dimensional  grid  generator  in  chapter  I.  That 

is,  a  weight  function  is  generated  from  the  available  data,  and  then  a  new  mesh 

(giving  the  x^  .'s)  is  generated  from  a  small  number  of  iterations.   With  the  flux 
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veloclties  calculated  at  the  vertices,  the  velocities  may  then  be  calculated 
everywhere  on  the  mesh  using  the  bilinear  interpolation  formula  (IV, 7). 

Next  the  boundaries  across  which  the  fluxes  move  must  be  determined.  The 
cell  centered  quantities  will  be  fluxed  using  the  mesh  lines  as  boundaries. 
However,  for  vertex  centered  quantities,  a  mesh  must  be  invented.  This  mesh  uses 
midpoints  between  vertices  and  cell  centers  as  the  corners  of  the  region  to  be 
fluxed  as  illustrated  in  Fig.  AA  by  the  dotted  line.  These  geometries  allow 
fluxes  to  be  computed  cell  by  cell.  Two  fluxes  are  computed  for  each  quantity. 
The  first  flux  is  generated  using  upwind  differencing  while  the  second  flux  will 
use  interpolated  donor  cell  differencing.  These  methods  are  well  known  and  are 
first  and  second  order  accurate  respectively  [1B,11B].  The  two  fluxes  are  then 
used  to  calculate  a  single  flux  via  Zalesak's  two-dimensional  FCT  algorithm  [13A] 
generalized  to  a  nonuniform  mesh. 

Mass  and  the  momentum  are  fluxed  between  vertices  and  int°nal  onpr^"-  19 
exchanged  between  cells.  As  a  result  mass,  momentum  and  internal  energy  are 
conserved  in  the  remap  phase.  However,  total  energy  is  not.  To  conserve  total 
energy  an  extra  step  is  added. 

Since  mass  and  momentum  exchange  are  done  on  a  cell  by  cell  basis,  the 
kinetic  energy  dissipated  between  four  vertices  of  a  cell  should  be  added  to  the 
internal  energy  of  the  cell.   To  do  so  the  following  formulas  are  used 


n+ln  n.   2 

^li         ^ii 

AKE.  ,  =  ii_  -  1^ (IV. A6) 

-■■J   o/n+l^   N    oz-n^  n+l„   \ 
2(   m^j)    2^  '"ij   ™ij^ 


where  '^q^^  is  the  momentum  at  the  ij  vertex  at  time  n  and  "m .  .  is  the  vertex 
mass.   This  can  then  be  rewritten  as       -   ' 


-81- 
n+l/2„   ,„      „  9 


where  dq  and  dm  are  the  changes  in  momentum  and  mass  of  a  vertex.  However  dq  and 
dm  are  the  sum  of  momentum  and  mass  fluxes  over  four  cells.  By  isolating  the 
momentum  and  mass  exchanges  for  one  cell  and  suciTning  over  the  resulting  changes 
in  kinetic  energies  this  sum  is  subtracted  from  the  internal  energy  of  the  cell. 
Thus  total  energy  is  conserved  to  roundoff. 

Numerical  Experiment  (Rayleigh-Taylor  Problem) 

The  Rayleigh-Taylor  problem  is  a  difficult  problem  to  do  with  a  Lagrangean 
mesh  because  of  mesh  tangling  that  takes  place  in  the  early  stages  of  the 
calculation.  To  do  this  problem  calculations  are  usually  done  on  an  Eulerian 
mesh.  In  the  numerical  experiment  presented  here,  the  mesh  is  caused  to  move  in 
a  manner  which  increases  the  resolution  of  gradients  in  the  solution  but  without 
tangling.  For  this  purpose,  the  weight  function  Vp/p  is  used  in  the  adaptive 
mesh  generator  with  X  set  to  1.0.  The  scaling  of  the  weight  function  Is  such 
that  its  value  ranges  from  one  to  one-hundred. 

The  experiment  will  be  set  up  on  the  computation  mesh  given  in  Figs.  A5  and 
46.  The  physical  size  of  the  computation  region  is  one  by  three  units  and  the 
density  jump  goes  from  one  to  two.  The  internal  energy  per  unit  mass  is  equal  to 
10  over  the  entire  domain  and  is  large  enough  that  the  sound  speed  is  very  large 
compared  with  fluid  velocities  as  in  incompressible  flow.  The  dark  line  across 
the  center  of  the  computation  mesh  is  formed  by  a  large  number  of  tracer 
particles  that  move  with  the  interface  between  the  two  fluids.  The  velocity  is 
initially  set  to  zero  and  the  gravitational  acceleration  set  to  10  in  the 
downward  (-y)  direction.   The  boundary  conditions  are  such  that  in  the  Lagrangean 
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phase  no  acceleration  is  allowed  normal  to  the  boundary  while  in  the  remap  phase 
no  fluxes  are  allowed  out  of  the  region  and  normal  velocity  components  are  set  to 
zero.   The  initial  perturbation  of  the  interface  is  large  (25%  of  the  width  of 
the  box). 

Figures  47,  48,  and  49  below  show  the  evolution  of  the  instability  (mode  1) 
after  a  short  time.  Figure  47  is  the  adaptive  computation  mesh  and  Fig.  48  is 
the  corresponding  density  contours.  Figure  49  is  an  Eulerian  calculation  using 
the  mesh  in  Fig.  45.  One  apparent  difference  between  the  solutions  is  that  the 
adaptive  mesh  solution  is  more  diffuse  than  the  fixed  mesh  solution.  The 
enlargements  of  the  interface  show  more  clearly  the  difference  in  interface 
thicknesses. 

As  the  solution  evolves  further  into  a  nonlinear  regime  the  solutions  look 
more  and  more  different.  Figures  50  thru  53  show  the  instability  at  a  late  time. 
Figure  50  is  the  adaptive  mesh,  Fig.  51  is  the  adaptive  solution,  Fig.  52  is  the 
solution  on  an  Eulerian  mesh  and  Fig.  56  is  the  solution  on  a  refined  Eulerian 
mesh  (40  x  120).  Figure  53  allows  one  to  compare  accuracy  of  solutions  50  thru 
52.  Two  resolved  features  are  especially  prominent.  The  first  feature  that  is 
resolved  by  the  adaptive  mesh  is  the  Helmholtz  rollup  near  the  bottom  of  the 
solution.  Notice  that  in  the  adaptive  solution  the  contours  of  the  densitv  form 
and  inverted  mushroom  while  the  Eulerian  calculation  contours  are  more  diffuse. 
Notice  also  at  the  point  of  symmetry  a  Kelvin-Helmholtz  type  instability  forms  in 
the  adaptive  case  while  it  is  prevented  from  forming  by  the  increased  numerical 
diffusion  in  the  fixed  Eulerian  case.  Even  on  the  finer  Eulerian  mesh  the 
instability  isn't  prominent.  .     , 

As  can  be  seen  from  the  numbers  of  cycles  used  (NCYC),  the  adaptive  method 
runs  slower  by  slightly  less  than  a  factor  of  two  than  the  Eulerian  version. 
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Increaslng  the  adaptlvity  by  increasing  the  volume  weighting  functional 
contribution  would  cause  finer  resolution  around  gradients  and  an  Increase  in  the 
number  of  cycles  per  calculation  would  result  because  of  stability  conditions  on 
the  timestep.  However,  if  a  uniform  mesh  were  refined,  the  amount  of  work  would 
increase  quadratically  from  the  increased  number  of  mesh  points  and  linearly  from 
the  increased  number  of  timesteps  required  by  stability  conditions.  Thus  uniform 
refinement  of  a  mesh  to  increase  resolution  of  gradients  by  a  factor  of  two  would 
increase  the  work  by  a  factor  of  eight,  compared  with  the  factor  of  two  increase 
in  work  for  an  adaptive  mesh  observed  above  for  a  similar  increase  in  accuracy. 

Numerical  Experiment  (Mach  3.0  Wind  Tunnel  with  a  Step) 

The  wind  tunnel  with  a  step  is  a  problem  that  has  been  used  to  test  many 
finite  difference  algorithms  in  two  dimensions  [lOA,  HA,  12A] .  Here  we  will 
demonstrate  the  accuracy  that  can  be  achieved  using  the  adaptive  mesh  with 
standard  low  order  difference  equations.  The  physical  region  is  a  1  by  3  unit 
rectangle  with  a  step  0.6  units  from  the  inflow  boundary.  The  step  height  is  0.2 
units.  Figure  54  illustrates  the  region  of  interest.  The  flow  will  be  initiated 
impulsively  meaning  that  the  flow  field  everywhere  will  be  set  to  3  units  per 
second.  Also  the  thermodynamic  variables  are  uniform  initially.  The  density  p 
is  set  to  l.A,  I  the  internal  energy  is  set  to  1.786  and  y  the  ratio  of  the 
specific  heats  of  the  fluid  is  also  1.4.  The  initial  pressure  and  sound  speed 
are  both  unity.  Thus  a  mach  number  of  3  is  arrived  at.  Finally  the  inflow  is 
set  to  maintain  the  above  quantities  at  the  left  boundary. 

One  of  the  first  to  test  this  problem  was  Emery  [lOA].  His  purpose  was  to 
compare  four  finite  difference  schemes.  However,  he  only  discussed  the  problem 
briefly  and  displayed  solutions  that  had  rather  poor  structure.  B.  van  Leer 
[HA]  also  used  this  problem.   His  solutions  demonstrated  the  merits  of  using  a 
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higher  order  version  of  Godunov's  method  [IB]  (specifically  van  Leer's  method  is 
second  order  accurate).  Woodward  and  Coella  [12A]  have  extended  van  Leer's  ideas 
to  a  third  order  method  and  have  gotten  solutions  with  superior  structure  [16A]. 
The  solutions  they  have  obtained  are  a  benchmark  for  this  problem  and  with  their 
permission  I  have  included  their  solutions  in  this  thesis  in  Appendix  C.  Woodward 
and  Coella  plan  to  do  a  review  paper  incorporating  this  problem  and  comparing  ten 
schemes.  However,  the  boundary  conditions  for  the  review  paper  at  the  step 
corner  have  been  changed  in  a  manner  that  radically  changes  the  nature  of  the 
solution  globally  from  the  results  presented  here. 

The  finite  difference  calculation  presented  in  this  paper  will  start  with 
the  zoning  in  Fig.  55.  The  adaptive  calculations  will  have  a  different  mesh  at 
later  times.   The  mesh  in  Fig.  55  has  40  x  120  =  4800  cells. 

Before  describing  some  calculations  further  attention  must  be  paid  to  the 
boundary  conditions.  The  outflow  boundaries  for  this  problem  are  Interesting 
because  of  the  diverse  opinions  on  how  to  deal  with  them.  Roache  [HB]  has 
argued  that  the  boundary  conditions  at  the  outflow  region  influence  the  upstream 
solution  even  when  the  flow  is  supersonic  (which  It  always  is  In  this  problem). 
However,  our  computations  seem  to  show  that  the  outflow  region  doesn't  influence 
the  computation.  (Among  the  boundary  conditions  which  were  tried  were  normal 
derivatives  of  all  variables  set  to  zero  (Neumann)  and  linear  extrapolation  of 
all  quantities.) 

Another  boundary  condition  that  must  be  considered  is  how  to  treat  the  two 
corners  of  the  step.  As  was  mentioned  previously,  the  solution  is  quite 
sensitive  to  these  boundary  conditions.  For  these  calculations  the  corners  were 
rounded  using  the  following  prescription.  A  line  is  drawn  between  two  adjacent 
boundary  points  and  as  shown  in  Fig.  56  the  fluid  is  allowed  to  flow  in  a 
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direction  tangent  to  that  line.   It  Is  this  boundary  condition  that  makes  it 
difficult  to  compare  this  or  any  other  solution  discussed  directly. 

The  following  figures  show  the  evolution  of  the  problem  through  a  time  t  of 
four  units.  For  these  calculations  the  gradient  length  of  the  pressure  |Vp|/p 
was  used  as  a  weight  function.  It  was  found  advantageous  to  scale  the  weight 
function  linearly  with  the  distance  from  the  corner  of  the  step.  Additional 
resolution  is  not  needed  at  the  corner  because  the  mesh  bunches  there  naturally. 
This  is  because  of  the  same  electrostatic  analogy  as  in  the  third  chapter  with 
the  aircraft  fuselage.  The  computation  mesh  and  pressure  contours  are  shown  in 
Figs.  57-72.  They  should  be  compared  with  Woodward's  solutions  in  Appendix  C.  We 
first  discuss  the  evolution  of  the  problem  in  time. 

At  time  0.5  a  shock  has  formed  in  front  of  the  corner  and  is  expanding 
towards  the  wall  opposite  the  step  (Fig.  58).  The  region  spoked  by  the  contours 
is  called  a  rarefaction  region  and  the  flow  there  is  transonic.  Notice  how  the 
mesh  (Fig.  57)  tries  to  align  itself  parallel  to  the  shock  to  resolve  the 
gradients.  This  shows  the  behavior  of  the  mesh  is  locally  one-dimensional  in 
n;iture.  This  characteristic  (see  [9A])  allows  us  to  use  the  knowledge  gained 
from  one-dimensional  problems  in  predicting  the  behavior  of  an  adaptive  mesh  in 
many  dimensions. 

At  time  1.0  the  shock  is  attached  to  the  wall  and  extends  toward  the  outflow 
boundary  (Fig.  60).  The  mesh  (Fig.  59)  is  uniformly  contracted  all  the  way  to 
the  outflow  region  because  the  weight  function  is  scaled  with  the  distance  from 
the  step. 

At  time  1.5  a  second  shock  reflection  enters  into  the  picture  (Fig.  62).  A 
slight  stem  or  mach  stem  is  forming  at  the  first  reflected  shock  point.  Notice 
that  the  mesh  (Fig.  61)  has  no  trouble  resolving  the  reflection  at  the  boundary. 
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Thls  is  because,  with  the  exception  of  the  corners  and  In  a  neighborhood  of  the 
step,  the  boundary  points  are  adjusted  to  keep  the  boundary  lines  orthogonal. 

At  time  2.0  a  Mach  stem  is  clearly  seen  at  the  first  reflection  location. 
The  intersection  point  of  the  three  shocks  is  called  a  triple  point  (Fig.  6A). 
Also  the  shock  strength  has  Increased  in  the  vicinity  of  the  outflow  region. 

At  time  2.5  the  Mach  stem  is  more  apparent  and  the  shock  strength  near  the 
outflow  region  is  greater  than  before  (Fig.  66).  The  mesh  does  not  tangle  at  the 
outflow  corner  because  the  smoothness  functional  keeps  points  apart  (Fig.  65). 

At  time  3.0  a  second  Mach  stem  has  formed  behind  the  step  and  the  third 
reflected  shock  becomes  apparent  (Fig.  68). 

At  time  3.5  (Fig.  70)  the  second  Mach  stem  lengthens  and  a  small  maxima  is 
seen  behind  the  triple  point. 

At  time  4.0  (Fig.  72)  the  most  complicated  structure  obtained  in  the 
transient  phase  is  seen.  Notice  how  the  first  Mach  stem  has  moved  forward  since 
it  was  formed. 

Figure  73  shows  a  solution  at  time  4.0  with  a  fixed  Eulerian  mesh.  Most  of 
the  structure,  with  the  exception  of  the  maxima  behind  the  triple  point,  is 
there.  However  the  entire  structure  is  more  diffuse.  As  one  looks  down  stream 
the  shock  strengths  are  much  weaker  than  the  corresponding  strengths  in  the 
adaptive  solution. 

There  is  not  a  great  deal  of  difference  between  the  number  of  computation 
cycles  used  by  the  adaptive  mesh  and  the  number  used  by  the  Eulerian  mesh.  The 
reason  for  this  is  that  the  fine  mesh  at  the  step  dominates  the  timestep  through 
the  stability  condition.  As  a  result,  it  cost  very  little  more  to  compute  with 
an  adaptive  mesh  for  this  case.  However,  from  the  illustrations,  the  adaptivlty 
gives  much  cleaner  and  more  detailed  contours. 
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To  conclude  this  section,  two  points  should  be  made.  First,  the  adaptive 
code  has  done  very  well  in  resolving  the  structure  of  the  problem  (120  x  40  cells 
versus  Woodward's  high  order  calculation  with  150  x  50  cells).  In  comparing  the 
results  with  Woodward's  calculations,  only  minor  differences  in  the  detailed 
structure  of  the  solutions  appear.  Some  of  these  differences  may  in  fact  be  due 
to  differences  in  Woodward's  boundary  conditions  at  the  corner.  (As  mentioned 
before,  the  solution  is  quite  sensitive  to  this  boundary  condition.) 

The  second  point  that  should  be  emphasized  is  the  adaptive  mesh's 
performance.  The  mesh,  as  seen  in  the  previous  Illustrations,  has  followed  the 
structure  of  the  problem  well.  What  is  particularly  encouraging  is  that  the  mesh 
resolved  the  shocks  regardless  of  their  orientation  without  mesh  tangling. 

Conclusions 

In  this  thesis  we  have  developed  a  comprehensive  and  flexible  set  of 
techniques  for  generating  adaptive  meshes.  The  formulation  prevents  mesh 
tangling  or  folding  by  building  in  mesh  smoothness.  Meshes  can  be  reliably 
generated  in  one,  two,  and  three  dimensions  in  a  manner  that  allows  one  to 
control  the  mesh  in  a  predictable  way.  These  techniques  are  also  advantageous 
because  they  are  not  directly  coupled  to  a  particular  finite  difference  scheme. 

We  have  shown  analytically  and  numerically  that  when  these  mesh  generators 
are  used  adaptively  for  singular  problems,  the  increase  in  the  amount  of  work 
required  to  perform  a  calculation,  in  most  cases,  scales  linearly  with  the 
reciprocal  of  the  minimum  cell  siz.e.  For  a  fixed  mesh  to  achieve  comparable 
accuracy,  a  uniform  refinement  must  be  done.  For  an  n-dimensional  calculation 
the  work  scales  as  the  n+1  sr  power  of  the  reciprocal  of  the  length  of  the  side 
of  a  computation  cell  in  the  mesh.   The  scaling  is  such  that  the  adaptive  mesh 
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brings   accurate   three-dimensional   calculations  within   the   capabilities   of 
existing  computers. 

The  adaptive  mesh  method  may  be  viewed  as  an  alternative  way  of  achieving 
more  accuracy  in  problems  which  have  boundaries  which  prevent  differencing  by 
standard  methods.  Instead  of  introducing  higher  order  schemes  and  the 
accompanying  difficulties  with  consistent  boundary  conditions,  one  can  use  a 
lower  order  body  fitted  code  and  resolve  the  problem  by  adaptively  zoning.  With 
body  fitted  coordinates  the  boundary  conditions  are  simple  to  implement  since 
they  are  lower  order. 

It  is  also  interesting  to  speculate  which  directions  one  can  take  in 
developing  the  methods  further.  For  example,  the  question  may  come  up  whether  it 
would  be  worth  while  to  combine  a  higher  order  scheme  with  an  adaptive  mesh.  In 
principle  this  is  an  unbeatable  combination,  but  in  practice  great  care  will  have 
to  be  taken.  For  example,  a  standard  high  order  scheme  reaches  out  to  several 
points  on  a  computation  mesh  to  interpolate  to  high  order.  In  general,  an 
adaptive  computation  mesh  may  be  quite  skewed  causing  truncation  errors  which 
would  degrade  the  accuracy.  To  avoid  this  problem  compact  methods  must  be 
developed  for  nonrectilinear  meshes  to  avoid  reaching  out  a  great  distance  on  a 
computation  mesh. 


-89- 
Appendix  A 

Existence,  uniqueness,  and  smoothness  of  the  Mesh  Generation  Equation 

Although  (1.15)  Is  computationally  convenient,  using  x  as  a  dependent 
variable  simplifies  the  proof  of  existence  and  uniqueness.  Using  (1.12)  the 
Euler  equation  becomes: 


^    W^H;^    ^'-   ^^^^x)'  +  tw(x)/C,  -  0  (A.l) 


or 


^xx  "  9 (A. 2) 

2[(1  -  t%/   +  tw) 


such  that  * 

Ki^l)   =  1   C(x^,)  =  N. 

Equation  (A. 2)  gives  a  recursion  relation  from  which  the  smoothness  of  the 
solution  5 (x)  may  be  determined. 

To  prove  the  existence  of  (A. 2),  a  continuation  method  will  be  used.  The 
continuation  parameter  t  is  built  into  the  equation.  The  proof  will  be  divided 
into  three  parts: 

(Step  1)   A-priori  estimates  will  be  provided  to  bound  i     and  its  derivatives 
independently  of  t. 
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(Step  2)   (A. 2)  will  then  be  shown  to  have  a  Frechet  derivative  (first  variation) 

with  a  bounded  inverse. 
(Step  3)   Steps  1)  and  2)  will  be  used  to  show  that  (A. 2)  has  a  solution  for  all 

t  in  [0,1].  - 

We  state  the  main  theorem:  " 

Theorem  1:   Suppose  w(x)  is  continuously  dif ferentiable  and  bounded  away  from 
zero: 


w(x)  >    v^   >   0 


Then,  the  solution  to  (A. 2)  exists  and  is  unique  for  all  t  e  [0,1]- 
To  prove  the  theorem  will  require  some  machinery. 

Theorem  2   (Implicit  Function  Theorem)  Suppose  F(x,y)  is  a  mapping  of  some 

neighborhood  of  the  origin  of  the  product  of  two  Banach  spaces 

into  a  third:  T 


F(x,y):  X  ^  Y  *  Z  (A. 3) 

such  that 

F(0,0)  =  0  ■-:.■■  ■     (A. 4) 

: .  .  '-..'■ .  J ■    ■  '■ ..  ''■'■■■  •- ; -■■' -  ■'    ■  '.^   '-'  '■  •  ■- ■ "        '■      ' ' 

Also  suppose  F  has  a  Frechet  derivative  with  respect  to  x  at  (0,0).   Further  F 
may  be  written  as 

■  ■,.•■,   -^  J  :'i:  !.■■■■  i<'     :'!     '  •!  ^ '■   •■••:."•■;  ^  5  u  i  "  .• .  . 
F(x,y)  -  Ax  +  R(x,y)  ,  ,^^  .;  ;  jro,  .-.:.  .  ^^'^^ 

such  that 
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a)  R  is  continuous  in  x  and  y 

b)  lR(xj,y)  -  R(x2,y)I  =  0(e)|xj  -  X2I  for  Ix^l,  IX2I,  lyl,  <  e 

c)  A  has  a  bounded  inverse  A~^  defined  on  all  z.  Then,  for  lyl  small, 
F(x,y)  "   0  has  a  unique  solution  x  «  x(y)  near  the  origin. 

In  using  Theorem  2  the  following  correspondences  will  be  made: 

F(x,y)  is  a  3  tuple  composed  of  the  differential  operator  and  the 
boundary  condition. 


F(x,y)  A  (  "^ !^LJ^ ,   5(xi)  -1   ,   Uxfj)  -  n)  (A. 6) 

=    2[(I  -  t)5/  +  tw] 


X  is  the  Banach  Space  C  [x^.Xj^]  space  of  twice  continuously  dlf ferentlable 


functions  on  (xj^,Xj>j)  with  norm: 


IIC(x)ll  =  max  Id  +  max   |E  |  +  max  \E,\  (A. 7) 

X  X      ^        X      ''^ 


Y  is  the  Banach  Space  (Interval)  [0,1]  with  norm: 

II  til  =  It!  (absolute  value)  (A. 8) 

Z  is  the  Banach  Space  C[xj^,x„],  the  space  of  cont.   functions  with  norm: 


It  will  be  shown  that  near  the  origin  the  Image  of  r  is  a  subset  of  C[xj^,Xj^] 
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IIUx)ll  =  max   1^  I 

X 

crossed  with  t^.   That  Is  Z  -  C[x^,Xj^]  x  t^  x  jl 
The  first  variation  is  now  calculated  formally: 

Let  5  satisfy  (A. 2)  and  the  boundary  conditions.   Let  h  »  5  +  ng  where  n  Is  a 

2 
parameter  and  g  is  an  arbitrary  C  [xj,Xjg]  function.   Then  the  first  variation  can 

be  written. 


3F[h]  I     ^^'^    \ 
9n   |n=0   3n  |n=0 


where  t'  is  a  three  component  vector 


w^t(C  +  ng)x 

V,  =  (^  +  ^8)xx  +  o . 

1  2[(1  -  t)(C  +  ng)^-^  +  wt] 


V2  =  (C  +  ng)lx^  -  1 


V3  =  a  +  ng)L  -  1  . 


Then  the  components  of  the  first  variation  are 


3 VI  w^tg^  3w^t  Cx^  gx 


9n  ln=0    ^'^   2[(1  -  tX^^  +  wt]   2[(l  -  t)5x^  +  wt]2 


3  V2 
3n  |n=0 


(A. 9) 


g(xi)  , 


-93- 


3  vo 

-5 I   ^  =  gCXw^ 


Proof  of  Theorem  1 

(Step  1)   A-priori  estimates  will  now  be  found. 
Lemma  1   Suppose  C(x)  is  a  solution  of  (A. 2),  then  5   does  not  take  the  value 
zero  on  [x^ ,Xv] • 
Proof   Suppose  that  C^C^)  "  ^  ^or  Xj^  <  x  <  Xj^.   Then  using  (A. 2) 


^  tvJ,  „  dx 

tx=/   "L^ .  (A. 10) 

5  2[(1  -  tK,"*  +  wt] 


For  t  >  0  and  in  a  neighborhood  of  x  such  that  215^1  (1  -  t)  <  wt 


l^x'  ^  /   lUl  IwJ  dx 


X  '   '   X  ' 


Letting  Iw^^/wl  =  f(x)  >    0 

X 

f(x)  k^i  <  f(x)  /  k^i  f(x)  dx 


or 


R'  -  fR  <  0   where    R  =  /  f  ICx'  dx 

X 

-/_  f  dx 
letting  g  =  e   '^      >  0 


^   (gR)  <  0 
dx 


(A. 11) 
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Integrating  from  x: 

gR  <  0 

Which  implies  in  a  neighborhood  of  x 
R  E  0 

Finally  implying  in  a  neighborhood  of  x 


IC,I  E  0   . 


Using  (A. 10)  ^x  =  0  on  the  entire  interval.   But  this  implies  E,    =  constant  which 
is  a  contradiction.   For  t  =  0  the  same  contradictions  results  from  a  direct 
calculation 
Lemma  2   Suppose  ? (x)  is  a  solution  of  (A. 2),  then  ^^  >  0. 
Proof  By  the  mean  value  theorem  there  exists  x^  such  that 


C^(Xq)  =  J! L  >  0   .  (A. 12) 


X 


N 


Since  ^   cannot  cross  zero  E,  ^   >  0  everywhere. 
Lemma  3  ^  is  bounded  independently  of  t.  In  fact  1  <  5  <  N. 

Proof   Suppose  there  exists  x^  such  that  5(Xq)  >  N,  Xj^  <  x^  <  Xj^  ,  then  by  the 
mean  value  theorem  there  exist  x,  such  that 


5(Xo)  -  N 


''o  "  ^H 
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A  similar  proof  holds  for  an  x„  such  that  £(x  )  <  1. 

o  ^  ^  o 


Lemma  4  ^^  ^^  bounded  Independently  of  t. 


Proof   By  the  mean  value  theorem  there  exists  a  point  x  such  that 


r  /   ^     N  -  1 

^x^'^o)  =  ^   ■  ^   "  Cl 

x^    Xj_ 


Then  using  (A. 2)  and  assuming  t  >  0, 


•,■(■•:  '  . ■ 

^  C  w  t 

^x  x"- 


x^   2[(1  -  tK^^  +   wt] 


(A. 13) 


X 

<  max  I  .;_  I  J    5   dx  +  ci 
V     2w        '^       •■■ 
x„ 


w 
<  max  \   —    \    (N  -  1)  +  c, 
X    2w  ^ 

(If  t  =  0  we  obtain  cj^) 

Lemma  5   Suppose  5  is  a  solution  of  (A. 2)  then  ^         is  bounded  independently  of  t 
Proof 
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Iw^^tl  ,. 


2l[(l  -  tK„-*  +  wt)|  -.. 


(A.U) 


N  -  1 
X  '  2w  '  ~x"  2w   ^"'    "   ^N  "  '^l 


<    max  — - —  (max-— —  (N  -  1)  + 


Lemma  6  The  solution  of  (A. 2)  Is  unique. 

Proof   Suppose  there  are  two  solutions  5  ,ti  of  (A. 2).   Then  letting  v  =  C  -  n 


t'^x^x  tw^n^ 


V   =  - ■ ■ 

2[(1  -  tK^3  +  wt]   2[(1  -  t)nx^  +  wt] 


[2t^w^  -  2tw^(l  -  xK/i^(Cx  +  nx))Vx 
2[(1  -  t)E,J   +  .L][(l  -  t)n^3  +  vt] 


(A. 15) 


=  f(x)v^ 


Since  v(x,)  =  v(xm)  =  0  by  Rolle's  theorem  there  exists  x,  x,  <  x  <  Xx,  such  that 
v^(5)  =  0. 

X  ,:'■■-,:- 

Vjj(x)  =  /  f(x)v(x)  dx  (A. 16) 


Using  the  same  technique  as  in  Lemma  1  v^(x)  =  0.   Then  using  the  boundary 
conditions  v  =  0. 
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(Step  2)  The  Frechet  derivative  (A. 9)  Is  meaningful  since  i^   is  bounded  away  and 
above  zero.   Since  the  boundary  conditions  for  (A. 2)  are  linear,  the 
boundary  conditions  of  (A. 9)  are 

gCx^)  =  0  -  g(x^)  (A. 17) 

(A. 9)  is  abbreviated. 

^[g]  =  gxx  "  P('^)8x   •  (A. 18) 


It  is  now  shown  that  the  Frechet  differential  operator  is  a  bounded  operator  from 

2 
C[x^,xj^]  to  C  [xj^.Xf^].   That  is  the  inverse  of  (A. 18)  with  boundary  conditions 

(A. 17)  is  a  bounded  operator  from  C[x^,Xjq]  to  C  [x^.Xj^].   To  show  this  we  can 

explicitly  write  down  the  solution  to 

L[g]  =  r(x) 

1  •   -  '  ■  V    .  i!  r. 

r 

with  zero  boundary  conditions. 


g(x)  =  Ci*l(x)  +  Ci.^j^'^^  "^  /   7^-^ (A. 19) 

w(5  ) 


/   (*ia)4'2'(x)  -  <l>2(5)<l'r(x)]  r(C)d5 
g'(x)  =  ci<t>i'(x)  +  C2't>2'(x)  +  /   7-^;^ (A. 20) 


g"(x)  =  ci4)]^"(x)  +  C24'2"(x)  +  r(x) 


(A. 21) 


+  / 
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'^l 


1^(0 


where 


^1 


''n    [^^(0(t.2(xN)   -  (^2a)*i(xN)]   ra)d5 

C2   =   -<fi(xi)   /  wa)[<^i(xiH2(''N)    -  'l>2(xi)'t>i(^N)] 

''I 


w(0    =  <t>i(0*2'(^)    -  <t'i'(^)<t'2(^) 


and 


/%(Od5 


<^l{x)   =  J       e  ds,        4i2(x)    =   1      . 

Using 

X  X 

I   /      g(x)r(x)   dxl   <    max    lr(x)|   /       lg(x)|    dx 

X 

x^  x^ 

and  noticing  everything   in   (A.19-A.21)    is  well   behaved 


max    (Igl    +    Ig'l    +    lg"l)  <    M,    max    |r|       .  ^  (A. 22) 

X  X 


Thus  by  Theorem  2  if  a  solution  (A. 2)  exists  for  t  »  t^,  then  in  a  neighborhood 
of  t  a  solution  exists  as  well.  Alternatively,  the  solution  of  (A. 2)  exists  for 
open  sets  on  [0,1]. 
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(Step  3)   We  next  show  that  If  <C(x,t^)>  is  a  sequence  of  solutions  such  that 

t^  •*■    t.   Then  there  exists  a  soln  5(x,t).   This  Implies  that  the  set 

S  =  [t:5(x,t)  is  a  soln  of  (A. 2)]  Is  closed. 
But,  there  are  only  two  subsets  of  [0,0]  that  are  open  and  closed.   They  are  (fi 
(null  set)  and  [0,1].   Since  a  solution  exists  for  t  =  0,1  then  a  solution  exists 
for  all  t  e  [0,1] 

Now  consider  a  sequence  <E.(x,t^)>   such  that  tj^  +  t  e  [0,1]. 

From  (step  1)  C,  ^'  are  equicontinuous  since  5'  and  E."  are  uniformly  bounded. 
Since  for  t  =  0  a  soln  exists  we  may  assume  t  >  0.  Then  all  but  a  finite  number 
of  the  t^   say  n  >  N  are  less  than  t/2.   Using  (A. 2) 

l"xx('^2)  -  "xx^'^l)! 

w^(x2)tj,Ux(x2)  w^(xi)tnUx(xi) 


2[(1  -  tj^)u^^(x2)  +  w(x2)  t^]    2[(1  -  t^)n^^(^i)    -  v(y.i)t^] 


tn 


2[(1  -  t^)u^3(x2)  +  w(x2)t„][(l  -  t^)uj(x^)   +  w(xi)t„] 


Iw^(x2)u(x2)[(l  -  tj,)u^3(x^)  +  w(x;^)t^]  -  w^(xi)u^(xi)[(l  -  tj,)u^3(^^)  +  "^(^l^^n^ 


for  n  >  N 


<  ? lu  rx^")n  fx^l  (I   -   t  ">n  /^ 


t^"min 


Wjj(x2)u^(x2)  (1  -  tj^)ux  (xj)  +  Wx^''2^"x^'^2^"(^l^^n 
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3/ 


-Wjj(xi)ux(xi)(l  -  tjj)ux-'(x2)  -  w^(xi)u3j(xi)w(x2)tj^l 


2  -3  -J 

< Iwjj(x2)ujj(x2)(l  -  tj^)Uj^^(xj)  -  Wj^(xj)Ujj(x;|,)(l  -  tjj)u^-'(x2)  I 

+  Iwj^(x2)u^(x2)w(x^)tj^  -  Wj^(xj)u^(xj^)w(x2)tjjl 


Now 

|w^Cx2)u../x2)(l  -  t^^Vi,.^(x^)  -  Wj.(x,Juj^(x.J(l  -  tj,)ux^(x2)l  (A. 23) 

<  Iw^(x2)  -  w^Cxj)!  !ux(x2)u^-'('x^)| 

+  Iu^(x2)  -  Ux^'^l)'  '"x^'^l^^x^^^l^' 

lu^(xj)  -  u^(x2)l  lwjj(xi)u^(xi)  I  lu^(x^)  -  u^(x2^)u^(x2)  +  u^^  (xj)! 

Since  w^  is  bounded  on  a  compact  interval  as  is  u^^ 

<  Iwx(x2)  -  w^(xi)|  Ml  +  Iux(x2)  -  Ux(xi)IM2 

Similarly 

Iwjj(x2)u^(x2)w(xi)tj^  -  Wjj(xj^)Ujj(xi)w(x2)tj^l  (A.2A) 

<  Iwj^(x2)  -  w^(xi)|M3  +  Iw(x2)  -  w^(xi)|M4  +  Iu^(x2)  -  Ux^^^l^'^S 

Finally 

'"xx('^2)  -  "xx^'^l)' 
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<  Mj'lu^Cxj)  -  "x^'^l^'  "^  M2'|w(x2)  -  wCx^)!  +  M3'|wj^(x2)  -  v^(x^)\ 

for  all  n  >  N 

Since  u  ,  w,  and  w  are  uniformly  cont.  on  the  compact  domain  the  set  of  u^^^'s 
for  n  >  N  are  equicontinuous.  But  the  rest  of  the  u^^'s  are  (n  <■  N)  are 
uniformly  continuous  on  compact  intervals  so  the  entire  set  ^^u^x^^-n^ '""^ ' '  ' ' '""^ 
is  equicontinuous.   Finally,  u(x,tj^),  u^(x,tj^),  and  u^^(x,tj^)  are  equicontinuous 

implying   subsequence   Ujj(x,t^  )   converging   uniformly   to   a   function   u(x,t) 

2 
uniformly  in  the  norm  of  C  [xj^.x^^].   u(x,t)  must  satisfy  the  differential 

equation. 

This  completes  the  proof. 
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Appendix  B 
Theorem  3   Suppose  v.",  for  all  meshes.  Is  bounded  everywhere  within  the  domain 
of  Interest  [0,1]  x  [0,T].   Then,  (I.A5)  to  (I.A7)  is  a  convergent 
scheme  provided 


c  +  max  II  v'^ll 

0  <  c^  <  [ ]  At  <  C2  (B.l) 

min  llh"ll 


max  h^" 

-J <  C3   .  (B.2) 

min  h." 


Proof   Let  the  error  be  written  as 


w."  =  u(x.",nAt)  -  *  j"   .  (R.3) 


A  difference  scheme  for  the  error  can  be  written. 


(c  -  v.^'^^^At  1  -  sgn(c  -  v.""^^^  _  2KAt 

J  _■----  ■J  ■  ■        —  .  -  — — —  -  — 


w-i+l    L  xi 2  ,    n+l/u   n+1  ■  V,  n+l^ 


n+l\A»  1  j_ /_    ..  n+l- 


„+l  r  -(^  -  Vi""')At  1  +  sgn(c  -  v^""M 


+  „   n+l  r  "^      :)    y"L  1  T  sgnvc  -  vj    ,    ^  2KAt 


V,  n+l  2  .  n+l/v,   n+l  .  v,  n+l\ 

h^  hj    (hj+i    +  hj    ) 


(B.A) 
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wj°"  1 


_^j        1   -  sgn   (c  - 


v,"-b 


At 


1   +  sgn    (c   - 


v-b 


^j+1 


n+1 


n+1 


(c   - 


vj"^b 


where 


+  ?fi! +1] 

V,  n+K        n+1  ■' 


=  Wj"  +  E"+l    (At)2 


1        X   >   0 


^'l      "   ^N      "   °      '         sgn(x)    =        0        X  =   0 


(B.5) 


-1      X   <   0 


and 


n+1    _ 


2  3^2  2         J  3^2  J  3t9x 


1   +   sgn(c   -  v.""^b  .,      h  ""^^  .2.         .. 

+  ^ (c  -  v^"+l)  4^  L4.  (X**.    (n+l)At) 


'J 


^t       3x2 


(B.6) 


1   -   sgn    (c   -   v.""^b 


n+1 


At  ^^2 
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[^ 


J+1 


n+l^  3^*  ,- 


3x- 


(X,  (n+l)At)  -  h  "+1  L*  (X,  (n+l)At)] 
■^    3x-^ 


+  tc 


3(hj+i"-^^  +  hj"+l)At 


The  starred  and  barred  quantities  are  intermediate  points  resulting  fron  the  use 
of  the  mean  value  theorem,  e"  ,  the  truncation  error  term,  using  (B.l),  (B.2), 
and  the  boundedness  of  the  derivatives  of  the  solution  of  (1.16),  is  bounded 
Independently  of  the  mesh  spacing. 

Let  j  be  the  index  of  w."  such  that  |w."|  =  II  w"ll  .  For  now,  assume  lw."|  *  0. 
Using  (B.A) 


iw. 


ni     |r    "j+l"  !     ^^   -   ^j"^^'^    ^    -    ^S"^^    -   ^j^ 


wy 


2KAt 


'j+1 


h^^l-Ch.^,"  .  hj") 


Wj_i"        -(c  -  Vj")At    1   +   sgn(c   -  Vj") 


"j 


2<At 


":" 


hj"(h3^i"  *  hj")  ) 


+   1  + 


ZicAt 


K  n+lu  n        V,        n 


''  1 '  - ""'; ' '"' )  (c  -  V,") 


(B.7) 


_^    ^t     [1   +   sgn(c   -   Vj")] 


^J^ 


(C     -    Vj")] 
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<  lwj"-M  +  |E"|  At2   . 


In  this  form,  It  is  easy  to  see  that  the  term  in  absolute  values  on  the  left  hand 
side  of  (B.7)  is  always  greater  than  one.   Hence 

llw"ll  <  llw""^ll  +IIE"ll(At)^   .  (B.8) 

By  induction 

n 

llw"ll    <    (At)2      I      IIE"iI 
j  =  l 

<  (At)^n  max  II  e"ii  (B.9) 

n 

<  (At)^   M  max  II  e""    <   At   T  max  II  e"ii 

n  n 

Since   It  E"ll    is   bounded 


|wj"l    <    kAt      .  (B.IO) 


Thus  the  scheme  is  convergent, 
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Appendix  C 


Woodward's  Third  Order  Solutions 
(Pressure  Contours) 


PR[SSU°[  PfiRpMJSL        S/ 

Dr5.e:L-C3    CGJRM  =  E.£C;!:     3Z  CONTOURS:  b.2z2i-Zl  TC 


:l*2\ 


0.5 


LD  CD  ^  C\J 


LD  CO  ^3  C^J  ZT"  LD  CD  ^ 

•  ••••••  " 

—  —  CsjCSjCNjCNjCNjr^ 


PRESSURE  PflRRMUSL      '9/32/80-  1      LMBSZPhC 

DT'4.E5:-23    COuRM-e.e01:     30  CONTOURS:  4.S64E-e!  TO    1.2:5E^B1 


t  =  1.0 
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FnLSSUR 


DT^5.C3E-E!3  CCJ^M-e.BC?:  30  CO\^GJBS:  5.  Hi!:-;;!  TO  l.Sc::'Cl 


t  =  1.5 


PRESSURE  PRRflMUSL   9/30/80-  1  IVS^mi 

D1-5.16L-23  CuunM-0,757:  30  CONTOURS:  4. 5Z2l-01  TO  1.251[*01 


t  =  2.0 
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D1-5.!9l-03    COURM  =  0.6t3:     3.^  C:\:3jRS:  4.67iL-Cl  10    \.\izi^2\ 


t  =  2.5 


PRESSURE  PflRflMUSL        9/3^/60-  1      LM352^nG 

D1-5.02E-03    C[!'5:,r0.e!4:     3Z  CCMCJRS:  4.e73[-ei  TO    l.lSSE^Bl 


t  =   3.0 
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E.e 


^  cxj  rr*  (_D  CD 


c^  z:-  (_o  CO 


CNJ 

=- 

OD 

CO 

!S 

C^ 

CXJ 

rj 

CSJ 

n~i 

PRESSURE  PhRR-^jSL 

N  •  705    T  ■  3  SC-5E'*2Z 
LCy,B3fl  SHOCK  PROBLEM.  E'J'.EPIRN  PflRR^'JSL.  C0URRNT-.8   S^XIEZ  GRID.  ^9/32/80 


3.5 


PRESSURE  PflRflMUSL        9/30/80-  1       LMB50nfll 

N  =     802         T  =  4.00260E+00 
LflMBDfl  SHOCK  PROBLEM.     EULERIR^  PflRflMUSL.     C0URflNT-.8      5i3X150  GRID.     9/30/80 


t  =  4.0 
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Flgure   Captions 

Figure  1     Steady  state  solution  of  equation  (1.16)  containing  a  boundary 

layer  near  x  =  1. 
Figures  2a, 2b  Analytic   solutions   and   derivatives   of  the  transient  part  of  the 

solution  of  equation  (1.16). 
Figure   3     Summary  of  the  one-dimensional  adaptive  algorithm. 
Figure   A     Infinitesimal  diagram 
Figure   5     A  graph  of  absolute  error  versus  Courant  number  for  the  uniform  and 

adaptive  schemes. 
Figure   6     Initial   conditions    for   a   typical   adaptive   one-dimensional 

calculation. 
Figure   7     The  adaptive  solution  at  t  =  .1  . 
Figure   8     The  weight  function  at  t  =  .1  . 
Figure   9     Initial  conditions  for  the  uniform  mesh. 
Figure   10    The  nonadaptive  solution  at  t  =  .1  . 
Figure   11    Absolute  error  is  plotted  as  a  function  of  the  grid  Reynolds  number 

(1.81)  for  the  uniform  and  adaptive  schemes. 
Figure   12     Illustration  of  the  parameter  space  in  two  dimensions. 
Figure   13    An  example   of   how   the   parameter  rectangle  can  be  mapped  onto  a 

nonsimply  connected  region  using  a  cut. 
Figure   14    A  mapping  of  the  parameter  space  rectangle  onto  the  physical   space 

rectangle. 

Figure   15    The  variance  over  the  mean  as  a  function  of  A  . 

2 
Figure   16     Contours  of  wJ   for  X^  =  0. 
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2 
Figure   17    Contours  of  wJ   for  X^  =  10. 


Figure   18    The  configuration  for  the  second  two-dimensional  experiment  mapping 

the  parameter  rectangle  onto  a  slit  disc. 
Figure   19    The   result   of   using   the   standard   smoothness   functional   for 

generating  a  mesh  on  a  slit  disc. 
Figure   20    The  mesh  collapses  when  periodic  boundary  conditions  are  used. 
Figure   21    A   randomly  displaced  polar  mesh  used  as  initial  conditions  for  the 

modified  smoothness  functional. 
Figure   22    The  result   after   20   iterations   using   the  modified   smoothness 

functional  on  the  mesh  in  figure  21. 
Figure   23    The  result  of  the  modified  weight  functional  with  a  radial  Gaussian 

weight  function. 
Figure   24    The  parameter  space  parallelepiped  in  three  dimensions. 
Figure   25    Mapping  of  the  parameter  space  to  a  cube  in  physical  space. 
Figure   26    Graph   of   the   values   of   the   three-dimensional   weiphtinc    ^nd 

smoothing  functional  values  as  functions  of  X^. 
Figures  27a, b  Cuts   in  physical  space  illustrating  how  the  mesh  is  deformed  using 

the  weight  function  (111.50). 
Figure   28    The  mapping  of  a  parameter  cube  onto  a  frustum  in  physical  space. 
Figures  29,30  Various  steps  in  generating  the  initial  mesh  for  the   orthogonality 

test. 

Figure   31    The  values   of   the   smoothness   and   orthogonality   integrals  are 

plotted  as  a  function  of  X,,. 

w 

Figure   32    Cuts  made  in  the  frustum  to  show  how  the  grid   tries   to   configure 
itself  in  a  way  that  minimizes  the  orthogonality  integral. 
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Flgure  33    The  upper  half  of  the  wing-fuselage  surface. 
Figures  34-37  A   sequence   of  sketches  showing  how  the  grid  Is  wrapped  around  the 

wing-fuselage  structure. 
Figures  38a, b  Cross  sections  of  the  Initial  mesh  configurations  for  the  winR   and 

fuselage. 
Figures  39a, b  The   same   cross   sections   after   50   smoothness   iterations  as  in 

figures  38a  and  b. 
Figure   40    The  support  of  the  weight  function  as  seen  from  the   front   of   the 

fuselage. 
Figures  41a, b  Cross   sections   of  the  mesh  after  50  iterations  with  the  weighting 

functional  having  a  nonzero  coefficient. 
Figure   42    A  schematic  of  the  piecewise  bilinear  representation  of  <f>  • 
Figure   43    Four  computation  cells  with  zero  velocity  at  the  boundary. 
Figure   44    An   illustration   of   the  mesh   and   fluxes   for  vertex   centered 

quantities. 
Figure   45    The  Eulerian  calculation  mesh  for  the  Rayleigh-Taylor  problem. 
Figure   46    The  initial  adaptive  grid  for  the  Rayleigh-Taylor  problem. 
Figure   47    The  adaptive  grid  at  an  early  time  (t  =  .5). 
Figure   48    Density  contours  for  the  adaptive  algorithm  at  t  =  .5  . 
Figure   49    Density  contours  for  the  Eulerian  algorithm  at  t  =  .5  . 
Figure   50    The  adaptive  grid  at  a  late  time  (t  =  1.8). 
Figure   51    Density  contours  for  the  adaptive  algorithm  at  t  =  1.8  . 
Figure   51    Density  contours  for  the  Eulerian  algorithm  at  t  =  1.8  . 
Figure   53    Density  contours  for  the  refined  Eulerian  scheme  at  t  =  1.8  . 
Figure   54    A  schematic  of  the  Mach  3.0  wind  tunnel  with  a  step. 
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Figure   55    The   initial   computation   mesh   generated   by   the   smoothness 

functional.  , 

Figure   56    Corner  flow  boundary  conditions. 
Figures  57-72  Meshes  and  pressure  contours  are  given  every  0.5  seconds  up  to   and 

including  t  =  A.O. 
Figure   73    The   solution   at  time  t  =  4.0  using  the  initial  mesh  fixed  for  the 

entire  calculation. 
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